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tic  e'ngi^nds  can  be  used  to  pump  heat  using  a  sound  wave  (refrigerator  or 
heat  pump)  or  pump  a  sound  wave  using  a  temperature  gradient  (prime  mover) .  The 
basic  arrangement  is  a  gas-filled  acoustic  resonator  with  appropriately,  positioned 
thenaoacoustic  elements.  Two  types  of  thermoacoustic  elements  are  used  in  these 
engines.  The  first  type  are  heat  exchanges  which  are  used  to  communicate  heat 
between  the  gas  and  external  heat ' reservoirs.  The  second  type  is  the  thermoacoustic 
engine  (lAE),  also  known  as  a  stack.  The  TAEs  are  sections  of  porous  media  that  sup¬ 
port  the  tec^erature  gradient,  that  transport  heat  on  the  acoustic  wave  between  the 
exchangers,  and  that  produces  or  Absorbs  acoustic  power.  Previous  theoretical  results 
in  thermoacoustics  have  been  developed  for  TAEs  with  circular  or  parallel  slit  pore 
geometries.  Ve  have  developed  a  general  linear  formulation  for  gas-filled  TAEs 
having  pores  of  arbitrary  cross-sectional  geometry.  This  analysis,  which  is  very 
helpful  in  designing  optimal  engines,  indicates  the  parallel  slit  pore  geometry 
optimizes  heat  and  work  flow.  Included  are  an  introductory  section  where  fundamen¬ 
tals  of  chermoacoustics  are  briefly  discussed,  a  summary  of  our  theoretical  analy¬ 
sis,  acoustical  measurements  of  an  air-filled  thermoacoustic  prime  mover,  and  numet- 
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ANNUAL  REPORT  ON 

THEORETICAL  AND  EXPERIMENTAL  STUDY  OF 
THERMOACOUSTIC  ENGINES 

Richard  Raspet,  Hemy  E.  Bass,  and  W.  Pat  Amott 
Physical  Acoustics  Research  Laboratory- 
University  of  Mississippi 
University,  Mississippi,  38677 

ABSTRACT 

Thermoacoustic  engines  can  be  used  to  pump  heat  using  a  sound  wave  (refrigerator  or  heat 
pump)  or  pump  a  sound  wave  using  a  temperature  gradient  (prime  mover).  The  basic  arrangement 
is  a  gas-filled  acoustic  resonator  with  appropriately  positioned  thermoacoustic  elements.  Two 
types  of  thermoacoustic  elements  are  used  in  these  engines.  The  first  type  are  heat  exchangers 
which  are  used  to  communicate  heat  between  the  gas  and  external  heat  reservoirs.  The  second  type 
is  the  thermoacoustic  engine  (TAE),  also  known  as  a  stack.  The  TAEs  are  sections  of  porous 
media  that  support  the  temperature  gradient,  that  transport  heat  on  the  acoustic  wave  between  the 
exchangers,  and  that  produce  or  absorbs  acoustic  power.  Previous  theoretical  results  in 
thermoacoustics  have  been  developed  for  TAEs  with  circular  or  parallel  slit  pore  geometries.  We 
have  developed  a  general  linear  fomulation  for  gas-filled  TAEs  having  pores  of  arbitrary  cross- 
sectiOnal  geometry.  This  analysis,  which  is  very  helpful  in  designing  optimal  engines,  indicates 
the  parallel  slit  pore  geometry  opdmizra  heat  and  wr^  flow.  Included  are  an  introductory  section 
where  fundamentals  of  thermoacoustics  are  briefly  discussed,  a  sumniary  of  our  theoretical 
analysis,  acoustical  measurements  of  an  air-Hlled  thermoacoustic  prime  mover,  and  numerical 
results  for  a  recently  constructed  helium-filled  prime  mover. 

INTRODUCTION 

Thermoacoustics  is  broadly  classified  as  the  interaction  of  heat  and  sound.  The  branch  of 
thermoacoustics  we  consider  is  heat  driven  oscillations  of  gas  in  a  tube  and  thermoacoustic 
transport  of  heat  The  basic  arrangement  is  shown  in  Fig.  la  where  thermoacoustic  elements  of 


heat  exchangers  and  a  TAE  are  shown  in  the  gas-filled  driven  resonance  tube.  The  plane  wave 
mode  of  the  resonator  is  considered  to  be  dominant  and  pressure  and  particle  velocity  have  near 
standing-wave  phasing.  This  arrangement  could  be  used  to  deliver  acoustic  power  to  the  TAE  for 
heat  transport  from  cold  to  hot  (for  low  temperature  gradients  (Th  -  Tc)/d)  as  in  a  normal 
refrigerator.  Thermoacoustic  elements  are  sections  of  capillaiy-tube-type  porous  media  as  shown 
in  Fig.  lb.  Theory  for  the  system  is  built  from  the  model  for  sound  in  a  single  arbitrary -perimeter 
capillary  tube  as  shown  in  Fig.  Ic.  Some  pow  perimeter  shapes  which  have  been  considered  are 
shown  in  Fig.  2.  Single  tube  radii  are  usually  designed  to  equal  the  frequency-dependent  thermal 
boundary  layer  thickness  5t  for  optimal  performance.  By  removing  the  acoustic  driver  and 
supplying  a  sufficiently  high  temperature  gradient,  the  TAE  produces  acoustic  power  at  the 
resonant  frequency  of  the  system. 

Observations  of  heat-driven  acoustic  oscillations  date  back  to  at  least  the  eighteenth  century. 
Rayleigh  f  gave  a  qualitative  explanation  well-worth  quoting:  "In  almost  all  cases  where  heat  is 
communicated  to  a  body,  expansion  ensues,  and  this  expansion  may  be  made  to  do  mechanical 
work.  If  the  phases  of  the  forces  thus  operative  be  favorable,  a  vibration  may  be  maintained." 
Acoustic  oscillations  were  noted  to  frequently  occur  in  a  capillary  tube  filled  with  helium  vapor 
with  one  end  of  the  tube  at  approximately  2  K  and  the  other  at  room  temperature  and  a  qualitative 
explanation  similar  to  Rayleigh's  was  given.^  A  full,  linear,  theoretical  investigation  of  heat-driven 
acoustic  oscillations  was  performed  fint  by  N.  Rott^  and  was  explored  in  a  series  of  papers 
starting  in  1969.  Rott  has  reviewed  this  work.^ 

The  reciprocal  mode  of  operation,  which  uses  a  sound  wave  in  a  resonator  to  transport  heat 
from  cold  to  hot  as  in  a  refngerator,  has  also  been  of  recent  interest.  This  thermoacoustic 
streaming  has  its  analogy  in  acoustic  streaming,  which  is  the  transport  of  mass  by  an  acoustic 
wave.^  Meikli  and  Thomann^  found  experimental  verification  for  their  theory  of  thermoacoustic 
streaming  in  a  driven  resonance  tube.  Wheatley,  Swift,  Hofler,  Garrett  and  othen  have  developed 
the  notion  that  the  arrangement  shown  in  Fig.  la  can  be  viewed  as  a  thermbdynamic  heat 
engine.'^'*  Swift  has  expertly  reviewed  this  work.^  The  thermodynamie  heat  engine  point  of  view 
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enhances  the  understanding  of  thermoacoustics  and  is  very  helpful  in  evaluating  practical  devices. 
The  connection  between  Refs.  4, 6,  and  9  has  been  briefly  explored  by  Rott.'^  Other  references  to 
the  eariy  history  and  practidoners  of  thermoacousdcs  can  be  found  in  the  review  ardcles.'^’^ 

A  common  approach  for  the  theory  of  sound  in  porous  media  is  to  envision  the  medium  as  a 
collection  of  circular  capillary  tubes. The  generalization  to  capillary  tubes  of  arbitrary  geometry 
has  recently  been  explored.  The  equations  and  boundary  conditions  used  in  porous  media 
modeling  and  in  thermoacoustics  are  nearly  identical  (thermoacoustics  has  an  extra  term 
proportional  to  the  ambient  temperature  gradient  that  occurs  when  evaluating  the  time  rate  of 
change  of  the  entropy).  It  was  apparent  to  us  that  thermoacoustic  theory  should  be  cast  in  a 
sufficiently  simple  form  that  analysis  of  pore  geometries  other  than  circles'^  and  parallel  slits'^’^ 
would  be  readily  possible. 

In  particular,  we  have  considered  use  of  extruded  ceramic  monolithic  catalyst  supports  (for 
example,  the  ceramic  used  in  some  automobile  catalytic  convertors)  for  use  in  thermoacoustics  on 
account  of  their  low  thermal  conductivity  and  regular  squaxs-pore  geometry.  This  material  has 
main  poits  of  diameter  154  mm.  In  addition,  die  walls  of  the  main  pores  are  porous  as  well,  with 
typical  wall  pores  of  diameter  25  )im.  The  frequency-dependent  complex  propagation  constant  of 
sound  in  the  square  pore  ceramic  Was  measured  in  related  woik.^'*  The  specific  acoustic 
impedance  of  a  49  cm  long  piece  was  measured  both  before  and  after  the  wall  pores  were  filled. 
These  measurements  verified  our  theory  for  sound  propagation  in  porous  wall  porous  media. 
The  TAE  used  in  the  lecture  demonstration  was  a  sample  of  this  ceramic. 

We  have  extended  thermoacoustic  theory  to  include  the  acoustic  field  quantities  and  the 
second  order  dnergy  flow  for  arbitrary  perimeter  pores.*^  Heat  and  work  flow  were  compared  in 
the  short  stack  approximation  to  investigate  the  eflects  of  pore  geometry.  Once  the  acoustical 
properties  of  the  separate  thermoacoustic  elements  have  been  determined  the  elements  must  be 
connected  in  series  inside  of  a  resonator  as  shown  in  Fig.  la.  Numerical  integration  of  the 
acoustical  equations  is  used  to  compute  field  quantities  in  the  stack  since  in  general  a  temperature 
gradient  exists  from  one  end  to  the  other.^*^  The  physical  parameters  of  ambient  density. 


viscosity,  sound  speed,  thermal  conductivity,  etc;  are  temperature  dependent  and  thus  depend  on 
location  within  the  TAE.  Specific  acoustic  impedance  and  pressure  translation  theorems  were 
developed to  compute  all  acoustical  field  quantities  and  energy  flow  at  each  point  in  the 
resonance  tube  shown  in  Fig.  la. 


THERMOACOUSTICS  IN  THE  LAGRANGUN 
AND  EULARUN  FLUID  FRAMES 


Figure  3a  and  3b  depict  an  inviscid  gas  parcel  oscillating  between  solid  paiallel  plates  at  the 
extremes  of  moticMi  in  the  Lagrangian  frame.  Heat  transfer  can  usually  be  neglected  in  normal  low 
frequency  acoustics.  However,  for  parcels  near  solid  surfaces  heat  exchange  is  likely  to  occur. 
Standing  wave  phasing  is  assumed  so  fluid  parcel  displacement  and  pressure  are  in  phase.  The 
plates  have  a  temperature  gradient  along  them  in  the  direction  of  oscillation.  The  location  for  wall 
temperature  To  is  nominally  the  parcel  equilibrium  location.  The  magnitude  of  the  temperature 
gradient  determines  whether  heat  is  transferred  to  or  irom  the  gas  parcel  The  right  end  of  the  plate 
faces  a  pr^ure  antinode  and  the  left  end  a  pressure  node.  Consider  fint  the  heat  pump  operation 
in  Hg.  3a.  The  gas  parcel  is  rapidly  displaced  left  a  distance  d  from  equilibrium.  During 
displacernent  parcel  pressure  and  hence  temperature  diminishes  and  expansion  occurs.  At  the  left* 


jtost  location,  the  parcel  is  momentarily  still  so  the  possibility  of  (slow)  heat  transfer  from  the  wall 


xurs.  For  a  small  wall  temperature  gradient,  the  parcel  is  now  at  a  lower  temperature  than  the 
all  so  Heat  flows  from  the  wall  to  the  parcel.  After  transfer  of  heat  dQ'  from  the  wall  to  the 
arcel,  the  parcel  expands  to  the  size  shown  by  the  shaded  square  and  wodt  is  done  by  the  parcel . 
a  the  gas  to  the  left.  This  work  discourages  vibration  of  the  gas  to  the  left  since  that  gas  is  in  the 
tparuion  part  of  the  cycle.  The  parcel  then  is  displaced  to  the  right  a  distance  2d.  Heat  dQ  is 
ansferred  to  the  wall  from  the  parcel,  and  the  parcel  shrinks  to  the  size  shown  by  the  shaded 


juare.  The  parcel  absorbs  work  from  the  gas  to  the  right  which  is  in  the  compressional  part  of  the 
roustic  cycle.  Net  beat  is  transported  by  the  parcel  from  the  wall  at  temperature  Tq  -  AT  to  the 
^all  at  temperature  Tq  AT.  It  is  now  easy  to  see  that  heat  from  an  exchanger  at  the  left  can  be 


transported, up  the  temperature  gradient  to  a  heat  exchanger  on  the  right  and  that  acoustic  power  is 
absorbed  by  the  parcel. 

Consider  next  the  prime  mover  in  Fig.  3b.  The  gas  parcel  is  rapidly  displaced  left  a  distance 
d  from  equilibrium.  For  a  large  wall  temperature  gradient,  the  parcel  is  now  at  a  higher 
temperature  than  the  wall  so  heat  flows  from  the  parcel  to  the  wall.  After  transfer  of  heat  dQ’  from 
the  parcel  to  the  wall,  the  parcel  shrinks  to  the  size  shown  by  the  shaded  square  and  work  is 
absorbed  by  the  parcel.  This  work  encourages  vibration  of  the  gas  to  the  left  since  that  gas  is  in  the 
expansion  part  of  the  cycle.  The  parcel  then  is  displaced  to  the  right  a  distance  2d.  Heat  dQ  is 
transferred  to  the  parcel  from  the  wall,  and  the  parcel  enlarges  to  the  size  shown  by  the  shaded 
square.  The  parcel  delivers  work  to  the  gas  on  the  right  which  is  in  the  compressional  part  of  the 
acoustic  cycle,  thus  encouraging  vibration.  Net  heat  is  transported  by  the  parcel  from  the  wall  at 
temperature  Tq  +  AT  to  the  wall  at  temperature  Tq  -  AT.  It  is.  now  easy  to  see  that  heat  from  an 
.  exchanger  at  the  nght  can  be  transported  down  the  temperature  gradient  to  a  heat  exchanger  on  the 
left  and  that  acoustic  power  is  produced  by  the  parceL 

Hgnre  4  shows  the  temperature  disturbance  for  AT=Q  at  a  fixed  coordinate  system  (Eulanan 
frame)  with  y-axis  as  shown  in  Fig.  3b.  The  excess  (or  acoustic)  temperature  changes  due  to 
pressure  changes  of  the  gas  and  conduction  reheat  to  the  walls.  The  assumed  boundary  condition 
for  Y»±l  is  that  fluid  and  wall  temperature  are  the  same.  The  quantity  shown  is  T(y)/Ta(i  = 
F(Y,Xt)  «  1  -  corAI(-i)l^  Xj  Y  /  2]  /  cosh[{ri)^>'^  Xy  /  2],  where  Y  -  y/(2a)  is  a  normalized 
coordinate,  Xt  ■  2a  (poo)Cp/K)f^  «  2^  a/Sx  is  a  normalized  number  proportional  to  the  ratio  of 
the  plate  spacing  2a  and  thermal  boundary  layer  thickness  5k:  Po>  Cp.  and  k  are  the  gas  ambient 
density,  isobaric  heat  capacity  per  unit  mass,  and  thermal  conductivity;  and  cu  is  the  radian 
frequency.*®  The  case  shown  has  Xy  -  3.2.  The  quantity  Tad  •  (Y  *  i)  Pi  /(P  Po  c^)  »  the 
acoustic  temperature  (driven  by  the  acoustic  pres'^ire  Pi)  that  would  occur  in  the  gas  if  no  heat 
transfer  occurred  as  in  normal  adiabatic  oscillation  of  sound  waves  (y  is  Cp/Cy,  p  is  the  coefficient 
of  thermal  expansion,  and  c  is  the  adiabatic  sound  speed).  The  magnitude  of  the  excess  gas 
temperature  nearly  reaches  the  adiabatic  sound  value  at  Y  a  0,  and  diminishes  to  zero  at  the  gas* 
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solid  boundary.  The  phase  difference  between  excess  temperature  and  pressure  occurs  because  as 
the  pressure  increases,  excess  temperature  increases,  aiid  particle  velocity  decreases  (on  account  ot 
standing  wave  phasing),  heat  is  transferred  to  the  wall,  cooling  the  gas.  The  excess  temperature  is 
a  maximum  at  a  time  between  maximum  particle  velocity  ai)d  maximum  pressure  at  a  given 
location.  This  time  phasing  between  motion  and  heating  is  a  common  feature  of  heat  engines.  The 
excess  temperature  for  non-zero  AT  is  modified  by  convective  transport  of  gas  and  by  gas 
viscosity.'^  Analysis  for  the  heat  transfeired  into  the  wall  when  AT  =  0  has  been  done.'^ 

Consider  the  single  pore  shown  in  Fig.  Ic.  Denote  the  first  order  acoustic  variables  as 
follows:  v(x,y,2)  =  v^(x,y,i)  +  V2(x,y,z)i'  is  particle  velocity,  Pi(z)  is  pressure,  ni(x,y,z)  is 
density,  and  Ti(x,y, 2)  is  temperature.  Ambient  quantities  are  indicated  by  subscript  0.  The 
lineaiued  fluid  equations  for  constant  frequency  oscillations  assuming  a  time  convention  erp(-ia)t) 

aiie^>9,l6 


-i  <fl  pO  vt(x,y,2)  i*  -  vt(x.y,2)  ,  (1) 

-i  0)  pi(x,y4:)  +  pQ^^)  Vx*vx(x,y,2)  =  0  ,  (2) 

pl(x,y,2)  a  -  po(2)  p  Tu'x,y,2)  +  iPi(z)  ,  (3) 

c*  ,  . 

-  i  CD  po(2)  Cp  Ti(x,y,2)  +  po(2)  Cp  Vt(x,y,2)  TOr  *  -i  co  p  To  Pi(2)  +  k  Ti(x,y,2)  ,  (4) 


where  the  transverse  gradient  and  L;»pla*.iari  operatews  are  defined  by  »  9/3x  H  t-  9/3y  P  and 
»  (32/dx2  +  92/8y2),  and  t|  is  gas  viscosity.  In  order,  these  equations  approximately  express  the 
2  component  of  the  equation  of  motion,  continuity  or  mass  conservation,  equations  of  state  for 
density,  and  heat  transfer.  Except  for  the  very  important  term  To.  ■  dTo/dx  in  Eq.  (4),  these  are 
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ths  equations  given  by  Zwikker  and  Kosten^l  in  their  solution  for  the  propagation  of  sound  in 
circular  pores.  Equation  (4)  expresses  that  the  temperature  at  a  fixed  position  changes  due  to 
motion  of  the  ambient  fluid,  due  to  compression  of  the  gas,  and  due  to  heat  conduction.  Boundary 
conditions  rse  v(x,y,z)=G  and  Ti(x,y,z)=0  at  ±e  pore  boundary. 

The  time  averaged  energy  flow  to  second  order  (subscript  2) 

H2(2)=^(z)  +  W2(z)  -  Qioss(z)  (5) 

where  time  ax'eraged  heat  flow  due  to  hydrodynanuc  transport  is 

^(z)  J(P0  ’^z(x.y.z)Ti*(x,y.z)  -  p  To  Vz(x,y.z)  Pi*(z))  dx  dy  ;  (6) 


the  heat  flow  by  conduction  down  a  temperature  gradient  is 

Qloss(z)  =  f*  Atts  "l*0t  +  (1  ~  fl)  Ares  ^tack  Tot  I  (7) 

and  the  time  averaged  work  flow  (or  power)  is 

W2(z)  ^  Jv,(x.y.z)  Pi*(z)  dxdy.  (8) 

Here  Ares  is  the  cioss'sectional  aien  of  the  resonance  tiibe  at  {»int  z,  A  is  the  cross-sectional  ^a  of 
a  single  pore,  Q  is  porosity,  as  <uid  Xstack  ^  tiic’  thermal  conductivity  of  the  gas  and  stack,  * 
indicates  complex  conjugation,  and  Ae  indicates  the  real  part  of  the  expression.  The  product 
(Q  Ares)  is  cross-sectional  open  area  of  the  tube  at  position  z. 
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Viscous  and  thennal  boundary  layer  thicknesses  are  given  by  5v  =  (2ri/copo)*^^  and  5^  = 
(2)C/a)poCp)l^‘.  We  introduce  a  dimensionless  "shear  wave  number"  X  =  R(poco/Ti)^'2  orX  =  2^'- 
R/5v,  where  R  is  a  characteristic  transverse  dimension  of  the  pore  in  Fig.  Ic,  and  a  dimensionless 
thermal  disturbance  number  Xf  =  R(pomcp/K)*^2  or  Xf  =  2^2  R/S^.  Use  of  the  Prandtl  number 
Npr  =  TiCp/K  gives  the  relation  Xt  =  X  por  definiteness,  take  R  to  be  twice  the  ratio  of  the 

transverse  pore  area  to  the  pore  perimeter  so  for  a  circular  or  square  pore,  R  is  the  pore  radius. 
Consider  the  newly  named  single  pore  transpon  function  F(x,y,X)  defined  by  the  following  partial 
differential  equation: 

F(x,y;X)  +  ^  F(x,y;X)  =  1  ,  (9) 

iX-* 


subject  to  the  boundary  condition  that  F(x,y,X)  =  0  at  the  (arbitrary  perimeter)  pore  wall  in  Fig.  Ic. 
In  anticipation  of  later  developmenc,  averages  over  the  pore  cross-section,  defined  for  example  by 
F(X)  ■  A*^  J  F(x,y;X)  dx  dy  where  A  is  the  area  of  the  pore  cross-section;  are  taken.  Acoustic 
particle  velocity  (1)  and  temperature  (4)  are  given  in  terms  of  pressure  and  F  by  the  relations 


vz(x,y,z) 


F(x,y;X)  dPi(z) 

itopo  ^ 


(10) 


Ti(x,y,2) 


v~  1 
po  P 


F(x,y;XT)  Pi(z)  - 


Tqi  .  F(x.y;XT)  ~  Npr  F(x,y;X)  dPt(z) 

podj'Z  1  -  Npr 


.(11) 


The  equation  for  pressure  is 


PO  dfF(X)dPi(z)>) 

F(X)d2[po  dz  J 


+  2a(X,XT)^^  +  k(X,XT)2  Pi(2)  =  0, 


(12) 


where 


1  -  N 


pr 


(13) 


and 


k(UT)2  °  ?  Pm  ‘'**■’■>1 


(14) 


In  the  absence  of  a  temperature  gradient  Tqz  =  0  so  ^  0.  The  complex  wavenumber  in  the 

pore  is  then  given  by  ±  k  which  is  the  usual  form  found  in  porous  media  modeling.  1 1 

SPECIFIC  ACOUSTIC  IMPEDANCE  AND 
PRESSURE  TRANSLATION  THEOREMS 

llie  specific  acoustic  impedance  of  an  acoustical  medium  is  equal  to  the  ratio  of  the  total 
acoustic  pressure  and  toted  particle  velocity.  The  acoitsdc  impedance  is  equal  to  the  ratio  of  the  total 
acoustic  pressure  and  the  total  volume  velocity.  For  porous  media  the  appropriate  boundary 
condition  at  a  surface  ate  continuity  of  pressure,  and  continuity  of  volume  velocity  or  equivalently 
continuity  of  acoustic  impedance,*®  Volume  velocity  is  At**  =»  Ares  Vz  where  £2  is  the 
porosity,  Vz  is  the  actual  particle  velocity  in  a  pore,  and  V2b  is  the  average  volume  velocity  per  unit 
area  of  pdrous  sample.  In  the  analysis  of  thetmoacoustic  elements  with  many  pores,  Vz  is  replaced 
with  Vzb  /  £2  m  all  of  the  single  pore  equations.  At  boundaries  between  thermbacoustic  elements, 
Pl(t)  and  Vj5(z)  or  equivalently  Z(z) «  P(z)  /  Vzb(2)  are  continuous,  where  Z(z)  is  the  specific 
acoustic  in^)edance. 

Rayleigh  developed  an  impedance  translation  theorem  for  homogeneous  fluid  layen.*'^  The 
impedance  translation  theorem  relates  the  specific  acoustic  impedance  at  one  side  of  a  layer  to  that 
at  the  other.  In  this  manner  one  may  apply  the  theorem  as  many  times  as  necessary  to  compute  the 
specific  acoustic  impedance  at  any  surface  in  the  layered  media.  This  translation  theorem  is 
applicable  to  heat  exchangers  and  resonattxr  sections,  but  is  not  applicable  to  the  stack  because  the 
physical  paiameten  such  as  density,  sound  speed,  ete,  depend  on  z  in  r  continuous  manner  mi 


account  of  the  temperature  gradient  Impedance  and  pressure  translation  theorems  which  take  into 
account  the  dependence  of  physical  parameters  on  position  were  derived  for  the  stack.  The 
expressions  are 


is  analogous  to  the  intrinsic  or  characteristic  impedance  of  a  porous  medium.  Expressions  (15)  and 
(16)  are  a  set  of  coupled  first  order  differential  equations  which  can  be  readily  solved  using 
numerical  techniques.  The  well-known  fourth  order  Runge-Kutta  algorithm  is  a  recommended 
numerical  method.  With  given  values  of  Pi(z)  and  2^z)  at  position  z  then  pressure  and  specific 
acoustic  impedance  are  determined  at  z-d  from  use  of  the  algorithm  Pifz-d) »  RK[Pi(z),Z(z)]  and 
Z(z<d)  «  RK(Z(z)]  where  RK  is  symbolic  notation  for  the  Runge-Kutta  algoridun. 

For  thermoacoustic  elements  not  having  temperature  gradients  (such  as  heat  exchangers  and 
sections  of  the  resonator),  Rayleigh’s  impedance  translation  theorem!^  can  be  written  for  porous 
mediaas 


2(z-d)  «  Zint 


Z(z)  coj(kd)  -  i  Zint  rin(kd) 
Zint  c<?^(kd)  -  i  Z(z)  rinfkd) 


(18) 


The  pressure  translation  theorem  is 


and 

W2(2)=:^  ReTiz)  .  (21) 

HEAT  AND  WORK  FLOW  IN  THE  SHORT  ST  XK 
APPROXIMATION  FOR  VARIOUS  STACK  GEOMETRIES 

The  short  stack  approximation  was  devised  by  Swift^  to  get  an  interpretablc  analytical 
expression  for  energy  flow  using  boundary  layer  theory.  For  irbitraiy  pore  shapes,  see  Ref.  16. 
Figure  la  shows  the  arrangement  for  the  short  stack  approximation.  Heat  exchangers  are  taken  u: 
be  of  negligible  thickness  and  thus  not  to  affect  near>8tanding  wave  phasing.  The  TAE  (or  stack) 
of  length  d  is  assumed  to  be  short  enough  that  the  empty  tube  standing  wave  is  marginally  affect^ 
The  temperature  difference  between  opposite  ends  of  the  TAE  is  assumed  to  be  much  less  than  the 
average  temperature  at  the  TAE  center  so  that  the  thermophysical  quantities  are  approximareiy 
constant  and  are  evaluated  at  the  average  temperature.  TAE  porosity  is 


With  a  rigid  tennination  at  z=0  in  Fig.  la.  specific  acoustic  impedance,  pressure,  particle 
velocity  amplitude,  particle  velocity,  and  particle  displacement  amplitude  and  particle  displacement 
at  z  are,  from  Eq.(  18)  and  Eq.  (19), 


Z(z)  =  -i  po  c  cot  koz  ; 


Pl(2)  -  PUCi)  coj  koz  ,  Vj(2)  =  ^f^ri/ikoz  ,  Vj(z)  =  ivzCz) , 

Qpoc 


(23a,b,c,d,e) 


where  the  wavenumber  in  the  empQ'  tube  is  ko,  the  n  phase  between  pressure  and  particle 
displacement  is  due  only  to  the  choice  of  the  coordinate  system  in  Fig.  la,  and  the  actual  pore 
particle  velocity  and  displacement  are  shown.  When  z  is  less  than  one  quarter  of  a  wavelength, 
Pl(z),  Vi(z),  and  ^x(z)  are  all  greater  than  zero.  Use  of  (22)  and  (23)  in  Eq.  (20)  gives  the  heat 
flow  at  the  hot  end  of  the  stack. 


Z  .  ,  QAres  Pl(z)  v|(z)  /m  [F*(Xt)/F{X)] 
Q2W  - - j - -  PTo 


PO  Cb  OAres  Vz(z)  /m(P*(Xr)  +  Npf  F(X)}  .s.  ,  Th  Tc 
2  (1  -  Npr2)  |F(X)|2  ^  d 


The  first  term  is  due  to  conversion  of  acoustic  power  to  heat  and  this  heat  flows  (in  the  TAE) 
towards  the  nearest  pressure  antihode.  The  second  term  is  heat  transported  on  account  of  the 
temperature  gradient  and  this  heat  flows  in  the  direction  t^tposite  to  the  positive  temperature 
gradient  irregardless  of  the  position  of  the  stack  in  the  standing  wave.  The  terms  preceding  the 
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temperature  gradient  are  kind  of  a  dynamical  coefficient  of  thermal  conduction.  The  TAE  acts  as  a 
refrigerator  when  the  first  term  is  larger  than  the  second.  Denote  by  Tqz  =  (Tc  -  Tnl'd  the 
temperature  gradient  across  the  stack,  and  denote  by  T  the  dimensionless  critical  temperature 
gradient  ratio  given  by 


r  3  To  c  tan  (kpz)  _  Cn  tan  (koz) 

””To  <oQ(y-1)  ^toapToc 

In  the  inviscid  approximadon  for  which  Nrjr  »  0  and  F(X) «  1, 

^(z)  « -  ^  Mi  0  To  f-(Xt)  (1  -  n  .  (26) 

Physically,  the  term  Im  F*(Xt)  is  a  measure  of  the  dynamical  thermal  interacdon  between  the  gas 
and  solid.  If  r<l  heat  is  transported  from  cold  to  hot  for  21co?’>  <tc.  According  to  Eq.  (26),  stacks 
made  of  pores  for  which  Im  F*(Xt)  is  a  large  value  will  result  in  the  greatest  heat  flow. 

Work  flew  is  given  generally  by  Eq.  (21).  No  woric  is  done  in  the  region  to  the  left  of  the 
stack  in  Fig.  la  because  in  this  region,  pressure  and  velocity  have  standing  wave  phasing.  To 
compute  the  work  done  in  the  stack,  use  is  made  of  the  impedance  transladon  <h?orem  to  get  the 
impedance  at  the  right  side  of  the  stacL^^  Denote  by  Vg  »  Ans  D  d  the  ambient  volume  of  gas  in 
the  TAE.  After  some  nunipdladon,  work  flow  to  first  order  in  kod  is 


WjW  tt-  D/mF-ftT)  - 

2  PO  c2  ^  1F(X)|2 


OA, 


res 


Pl(2)  vl(z)  P  (Th-Tc) 


/m{F*(XT)/F(X)} 


1  -N 


(27) 


pr 
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The  first  and  second  terms,  always  <  0,  are  dissipation  of  potential  and  kinetic  energy  per  unit  time 
due  to  thermal  and  viscous  diffusion  processes.  The  third  term,  which  is  >  0  when  the  hot  end 
faces  a  pressure  antinode,  is  the  acoustic  power  produced  on  account  of  the  temperature  gradient. 
When  the  third  term  is  larger  than  sum  of  the  first  two,  the  TAE  produces  net  acoustic  power. 
Woiic  flow  in  the  inviscid  approximadon  is 

W2(z)  Pi (0)2  ~  -  fm  F*(Xt)  cor2(koz)  (1  -  f)  .  (28) 

^  PO  Cp  . 

TAEs  made  of  pores  for  which  Im  F*  (Xj)  is  a  large  value  will  result  in  the  greatest  work  flow. 

Work  and  heat  flow  are  to  be  compared  for  the  various  pore  geometries  shown  in  Fig.  2a-2d. 
In  the  inviscid  short  stack  approximadon,  pores  with  a  la'ge  value  of  Im  F(X)*  will  have  the 
greatest  heat  and  work  flows  as  indicated  by  Eq.  (26)  and  Eq.  (28).  According  to  Fig.  5,  which 
shows  the  real  and  imaginary  parts  of  F(X)  for  the  various  pore  geometries,  the  parallel  plate 
geometry  has  the  largest  value  of  Im  F(K)*.  The  value  occurs  for  Xc  *  3.2  which  allows  one  to 
compute  the  optimal  operating  flequency  from  the  reladon  X;  >  (po  co  Cp  /  k)  R.  In  other  words, 
you  can  get  about  10%  more  heat  flow  and  work  flow  in  themoacousdes  by  choosing  to  make 
your  stack  from  parallel  plates  rather  than  the  other  pore  geometries.  The  funcdonal  form  of  F(X) 
for  the  various  pore  geometries  is  given  in  Ref.  16. 

DEMONSTRATION  AIR-FILLED  THERMOACOUSTIC  PRIME  MOVER 

The  air-filled  thermoacousdc  source  demonstrated  during  the  lecture  is  shown  schemadcally 
in  Fig.  6.  Heat  exchangers  wme  parallel  plates  of  copper  and  the  TAE  is  a  monolithic  catalyst 
support  extruded  ceramic.*2-i3  The  two-microphone-technique  impedance-tube**  at  the  bottom 
was  removed  for  the  demonstradon.  The  nominally  quarter-wavelength  resonator  was  capped  at 
the  top  with  a  rigid  plate  and  was  open  at  the  bottom.  Heat  tape  was  wrapped  around  the  hot  end 


(the  region  including  the  tube  at  top  and  the  first  heat  exchanger)  followed  by  heat  insulation. 
Water  was  circulated  around  a  jacket  surrounding  the  cold  heat  exchanger.  In  this  way,  i 
temperature  gradient  was  established  across  the  TAE,  For  sufficiently  high  temperature  gradients, 
(AT  =  176  K)  the  tube  produces  sound  at  116  Hz.  The  ambient  temperature  Tc  =  295  K.  The 
location  z=0  is  nominally  a  pressure  node  and  particle  velocity  antinode  so  the  specific  acoustic 
impedance  (equal  to  pressure  divided  by  particle  velocity  and  abbreviated  by  SAI)  is  a  relative 
minima  for  frequencies  in  the  vicinity  of  the  quarter  wavelength  resonance.  All  parts  were  made  of 
copper  except  the  TAE. 

SAI  measurements  were  made  as  a  function  of  the  temperature  gradient  across  the  stack.  The 
real  part  is  shown  in  Fig.  7  and  imaginary  in  Fig.  8.  Among  other  uses,  these  measurements  are 
helpful  for  evaluating  the  possibility  of  using  the  prime  mover  as  a  sound  source.'^  One 
interpretation  of  Fig,  7  is,  for  example,  that  the  plane  wave  reflection  coefficient  at  80  Hz  and  AT  = 
160  K  is  >  1  for  waves  incident  in  a  tube  of  the  same  diameter  as  the  prime  mover  but  in  the 
location  of  the  impedance  tube  in  Fig.  6.^  Also  shown  as  dashed  lines  in  Fig.  7  and  8  is  the 
expression  for  radiation  impedance  for  a  unflanged  tube.^l  The  expression  is  Zrad(co)  =  ~ 
pOC[(koR^)^  ^  i  ^'bere  ko  «  qVc.  To  determine  the  frequency  of  oscillation,  we  solve  for 

the  value  of  o>  such  that  the  complex  equation  Zrad(to)  *  Z(b})  measured  at  z=0.  The  operating 
point  AT  a  176  K  and  fres  «  116  Hz  of  the  prime  mover  demonstrated  is  shown  by  the  plus 
symbol  in  Figs.  7  and  8. 

To  validate  the  theoretical  model,  SAI,  denoted  by  Z(z«0;ci)c)  can  be  calculated  at  z=0  using 
(15)  and  (18).  Most  generally,  for  arbitrary  termination  impedance  Zend(ci>c)  at  z^O,  the  frequency 
of  operation  is  determined  by  setting  ZendCco^;) »  Z(z>«0;cOc).^  Here  cflt «  2nfres  -  i  Jcfres/Q  is  the 
complex  eigenfrequency  of  the  system  where  f(«$  is  the  resonant  frequency  and  Q  is  the  quality 
factor.  Notice  that  coc  is  a  function  of  the  tube  geometry  and  of  AT.  The  condition 
determines  AT  for  onset  of  acoustic  oscillation.  The  complex  eigenfrequency  has  recently  been 
investigated  for  a  helium  or  argon-filled  prime  mover  where  the  resonance  tube  was  sealed  at  z*0 
by  a  rigid  cap.23  a  word  of  caution*  For  prime  movers  above  onset,  or  for  strongly  driven 
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thennoacoustic  refrigerators,  the  temperature  distribution  from  hot  to  cold  is  not  determined  by  (7) 
alone^'^^,  and  use  caution  when  applying  (7)  to  determine  the  temperature  distribution  for  a  non- 
driven  tube  since  thermal  conductivities  depend  on  temperature.  The  presence  of  the  strong 
acoustic  wave  influences  its  thermal  surroundings^  by  heat  transport  and  in  this  sense  the 
thermoacoustic  oscillation  is  an  example  of  a  selfrinteracting  wave  process. 

With  the  impedance  tube  removed  in  Bg.  6,  measurements  were  made  of  the  sound  spectrum 
produced  by  the  prime  mover.  Figure  9  shows  the  spectrum  for  the  onset  temperature  AT  =  176  K 
(the  peak  at  425  Hz  was  not  produced  by  the  prime  mover)  and  for  AT  =  209  K.  The  ambient 
temperature  was  298  K.  A  Bruel  and  Kjaer  type  4147  1/2”  microphone  was  placed  at  z=0  in  Fig. 
6  for  these  measurements.  The  electrical  power  delivered  to  the  heat  tape  was  220  watts.  The 
value  AT  »  209  K  was  the  steady  state  equilibrium  temperature  gradient  of  the  system  where  heat 
supplied  by  the  upe  was  balanced  by  the  sum  of  the  acoustic  energy  radiated  away  and  dissipated 
in  the  resonator,  and  the  heat  deposited  at  the  cold  heat  exchanger  due  to  thermoacoustic  transport 
plus  normal  heat  conduction  dov/n  the  gradient  The  acoustic  power  radiated  at  1 16  Hz  and  AT  = 
209  K  was  esdmated  to  be  0.25  watts.  Some  of  the  applied  heat  energy  goes  into  the  acoustical 
energy  in  the  band  near  1 16  Hz,  in  the  higher  harmonics,  and  sturdy  DC  circulation  of  gas 
due  to  acoustic  streaming^^  that  occurs  out  of  the  prime  mover  at  the  center  and  returns  along  the 
walls.^  Evidence  of  acoustic  streaming  is  indicated  in  Bg.  9  where  the  generally  elevated 
background  level  is  due  to  flow  induced  noise.  A  flow  velocity  of  4  m/s  was  measured  for  the  gas 
exiting  the  central  part  of  the  prime  mover.  It  is  notewortiiy  that  the  eigenfrequencies  of  the  quarter 
wavelength  oscillator  are  given  by  fresC^m-i-l)  where  m>0,l,...  and  the  harmonics  of  the 

fundamental  due  to  nonlinear  processes^  are  given  by  fh  «  fns  n  where  n  m  2,3 .  For  example, 

the  peak  at  232  Hz  does  not  correspond  to  any  eigenfrequency  of  the  tube,  but  is  predicted  by  the 
analysis  of  strong  waves  in  open  tubes.^  Nonlinearities  in  closed  thermoacoustic  oscillators  have 
also  been  studied,  26 

Figure  10  shows  AT  «  209  K  sound  produced  by  the  prime  mover  with  the  microphone  in 
the  opening  and  5.7  cm  from  the  bottom  (inside  of  the  tube).  At  the  inside  point  the  sound 


pressure  level  is  145  dB,  1 16  Hz.  The  opening  cotresponds  nominally  to  a  pressure  node,  so  the 
acoustic  pressure  »  cos{xxJlV\  where  x  is  the  distance  from  the  top  and  L  =  72.7  cm  is  the  overall 
length  of  the  prime  mover.  Consequently  the  pressure  at  the  top  is  approximately  8  times  the  level 
at  5.7  cm  inside,  or  »  163  dB  at  1 16  Hz.  The  tube  radius  R  =  4.32  cm  is  much  less  than  the 
acoustic  wavelength  =  297  cm  at  1 16  Hz,  so  the  prime  mover  mouth  behaves  as  a  point  source  of 
spherically  expanding  waves.  Figure  1 1  sho'vs  the  sound  pressure  level  for  the  fundamental  at 
1 16  Hz  and  its  harmonics  5.7  cm  inside  the  prime  mover,  at  the  mouth,  and  at  distances  away  from 
the  mouth.  Efforts  to  increase  the  radiation  efficiency  of  the  prime  mover  in  a  given  direction 
would  result  in  more  radiated  acoustic  energy  and  would  need  to  be  accommodated  by  increasing 
the  temperature  gradient 

In  one  experiment  the  sound  production  by  the  prime  mover  was  suppressed  and  it  was 
super-heated  to  AT  »  285  K,  well  beyond  the  minimal  onset  AT  =»  176  K.  Figure  12  shows  the 
time  evolution  of  the  super-heated  prime  mover.  A  microphone  was  placfd  -  14  cm  from  the 
mouth  for  this  measurement  The  top  figure  shows  the  first  0.6  seconds.  Up  to  1.4  seconds  the 
amplitude  grows  exponentially  with  time  and  beyond  it  begins  to  level  off.  The  bottom  figure  is  a 
peak  detection  of  the  maximal  pressure  amplitude,  the  beginnings  of  which  is  shown  in  the  top 
figure.  Between  U  and  2  seconds  the  prime  mover  flutters  as  shown  on  the  bottom  figure  and  is 
apparent  to  an  observer.  During  this  time  the  built-up  heat  in  the  hot  end  is  used  up  by  sound 
production  and  thermoacoustic  heat  tianspon  down  the  temperature  gradient  After  2  seconds  the 
prime  mover  a  long  transition  to  a  steady-state  pressure  level  arxi  the  temperature  at  the  hot 
end  slowly  diminishes  to  the  steady  state  equilibrium  value  AT  •  209  K.  The  time  evolution  of  a 
non  superheated  prime  mover  has  been  studied  by  Muller,  et  aL^ 

NUMERICAL  RESULTS  FOR  A  HELIUM-FILLED  PRIME  MOVER 

Figure  13  is  a  summary  of  the  UM  thermoacoustic  engine.  Design  criteria  for  this  tube  was 
that  it  would  produce  sound  at  as  small  a  temperature  gradient  as  possible  given  the  constraints  on 
construction  of  parts.  After  construction  the  TAB  was  put  through  tests  to  see  if  onset  would  occur 


and  indeed  it  di±  The  tube  was  filled  mostly  with  Helium  lo  a  pressure  of  3  IcPa.  Appendix  A 
contains  the  programs  used  for  the  results  discussed  in  this  section. 

Figure  14  shows  the  computed  stability  curve  for  the  fundamental  frequency  near  308  Hz  and 
the  1st  harmonic  near  605  Hz.  Ambient  pressure  is  the  horizontal  axis.  Temperature  difference 
between  hot  and  cold  ends  is  on  the  vertical  axis.  For  temperatures  below  the  boundary  between 
sta*  le  and  unstable  no  gas  oscillations  occur.  Thermal  boundary  layer  thickness  5t  - 
(2ic/po<ocp)l^  <*  ari/2  where  Pq  is  the  ambient  pressure  in  the  tube.  For  fixed  tube  length 
as  in  the  UM  TAE,  the  boundary  layer  thickness  can  be  adjusted  by  changing  Pq.  For  low  Po,  the 
boundary  layer  is  much  thicker  than  the  pore  size  in  the  stack  and  heat  exchangers  so  gas  viscosity 
chokes  the  flow.  For  Pq  »  173  kPa  the  boundary  layer  thickness  is  optimal  for  thermoacoustic 
effects  and  the  lowest  onset  temperature  of  around  180  C  occurs.  One  factor  that  greatly 
contributes  to  miniiruzing  the  onset  temperature  is  the  Ideation  of  the  stack  and  heat  exchangers 
relative  to  the  lengths  of  the  hot  and  cold  ends.  When  the  elements  are  too  close  to  the  hot  end  gas 
particle  displacement  is  small  so  the  required  temperature  gradient  as  computed  from  ^  ;.  (25) 
becomes  prohibitively  large.  When  you  try  to  do  work  on  the  gas  it  responds  by  changing 
pressure  but  only  a  small  particle  displacement  occurs.  However  u  the  elements  move  too  close  to 
the  center  of  the  tube  where  a  pressure  node  occurs  the  work  done  approaches  zero.  When  you  try 
to  do  work  on  the  gas  at  this  point  in  the  standing  wave  the  gas  responds  by  undergoing  a  large 
displacement,  but  only  a  small  pressure  change  occurs.  Somewhere  betw  een  the  pressure  node  at 
the  center  and  the  particle  velocity  node  at  the  end  is  an  optinul  location  for  the  location  of  the 
elements.  The  optimal  location  is  closer  to  the  end  where  particle  velocity  is  less  so  that  losses  due 
to  gas  viscosity  are  minimal.  The  UM  TAE  was  designed  to  optimize  the  location  of  the  elements. 

Since  5t  •«  Tq®-*  tar  1/2  Pq^i/^  it  seems  odd  at  first  glance  that  die  minimum  of  the  first 
harmonic  would  be  at  a  higher  pressure  than  the  fundamental.  The  temperature  dependence  is  in 
part  due  to  ambient  density  and  in  part  to  thermal  conductivity.  Supposing  there  to  exist  an  optimal 
boundary  layer  thickness,  then  as  ci)  is  increased,  Pq  should  decrease  to  balance  the  equation. 
However,  there  is  no  single  optimal  boundary  layer  thickness  for  the  tubd.  Viscous  losses  can  be 


defeated  by  aiming  at  thin  boundary  layers.  The  stack  location  in  the  TAE  helps  to  minimize 
viscous  losses  for  the  fundamental,  bur  not  the  first  harmonic.  Additionally,  particle  displacement 
is  smaller  at  higher  frequencies  so  the  necessary  temperature  gradient  |Tozl  in  Eq.  (25)  to  make  F  > 
1  for  onset  increases.  (The  pressure  in  the  standing  wave  changes  from  maximum  to  minimum 
over  a  shorter  distance  as  the  frequency  increases.)  At  a  frequency  double  ihe  fundamental,  the 
ambient  pressure  must  diminish  by  a  factor  of  2*^2  to  get  the  same  optimal  boundary  layer 
thickness.  But  according  to  Eq.  (25)  the  average  temperature  (Th  +  Tc)/2  must  increase  by  a 
factor  of  2.  To  offset  the  temperature  increase,  the  ambient  pressure  must  increase  by  a  factor  of 
2^'^.  Thus  it  is  reasonable  that  the  ambient  pressure  must  increase  by  a  net  amouilt  approximately 
equal  to  79-^  to  achieve  the  bounclaty  layer  thickness  for  optintal  thermoacoustic  pumping  of  the 
wave.  In  order  to  experimentally  observe  both  modes,  the  TAE  should  be  heated  to  a  temperature 
in  the  unstable  region  of  the  1st  harmonic  with  the  Q  of  the  resonant  cavity  so  low  that  there  is  no 
oscillation  of  the  fundamental.  When  the  Q  is  increased,  both  modes  should  oscillate . 

Hguie  15a  is  the  1/Q  curve  for  the  fundamental  and  1st  harmonic  for  a  pressure  of  173  kPa. 
When  1/Q  a  0  the  TAE  is  at  the  perilous  boundary  between  stability  and  instability.  For  AT  above 
the  onset  value  154.6  K  for  the  fundamental,  the  tirrre  evolution  of  the  pressure  follows  the  form 
exp(Jifot/|Q|)  until  nonlinearity  becomes  apparent 27  The  exponential  growth  factor  for  the 
fundamental  has  a  max  imum  at  AT  «  800  K.  For  higher  temperatures  the  ambient  temper-ature  in 
the  hot  end  causes  the  ambient  souna  speed,  which  is  c  a  Tq^^,  to  increase.  As  a  rough  guide  for 
computing  the  resonant  fnequency  of  the  tube,  note  that  resonance  should  occur  for  2kocI< 
21co!iLh  ■  2rc  or  alternatively  fo  •  l/(2(Lc/Cc  Wch)]  where  Lc  and  Cc  are  the  lengths  and  sound 
speeds  of  the  cold  section,  etc.  Thus  as  Ch  increases  we  can  view  this  as  an  effective  shortening  of 
the  hot  end,  which  puis  the  elements  closer  to  a  velocity  node  with  concomitant  decrease  of 
effectiveness.  We  can  also  understand  the  increase  of  the  resonant  frequency  with  temperature 
shown  in  Fig.  15b  from  the  equation  for  fq.  Dispersion  effects  in  the  thermoacoustic  elements  are 
somewhat  apparent  in  this  figure  in  that  the  frequency  of  the  flrst  harmonic  is  not  twice  the 
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fundamentil.  This  difference  is  much  more  noticeable  in  the  original  numbers  from  which  these 
files  were  made. 

GENERAL  PRIME  MOVER  ANALYSIS  USING 
ENERGY  SOUNDARY  CONDITIONS 

We  have  analyzed  prime  movers  by  assuming  the  temperature  distribution  in  the  stack  from 
hot  to  cold  was  simply  of  the  form  To(z)  =  Tocoid  +  zfTohot  *  rocold)/d  where  d  is  stack  length. 
This  assumption  is  generally  invalid  for  two  reasons.  First,  even  for  temperature  gradients  below 
onset,  temperature  dependence  of  the  thermal  conductivities  of  the  stack  and  gas  cause  the  ambient 
temperature  distribution  to  differ  from  the  simple  form  above.  This  is  seen  in  Ref.  8.  Treating  the 
heat  flow  like  a  electrical  current,  the  applied  temperature  difference  at  both  ends  as  voltage,  and 
the  stack  and  gas  as  resistors,  net  resistance  is  the  parallel  combination  of  gas  and  stack.  The  gas 
and  stack  can  each  individually  be  considered  series  combinations  of  resistors  with  each  resistor 
h^ing  a  locally  isothermal  section  of  the  stack  or  gas.  Just  as  most  of  the  voltage  is  dropped  across 
the  largest  resistor  in  a  series  circuit,  the  largest  temperature  gradient  occun  in  the  stack  where  ute 
thennai  conductivity  is  lowest  For  gases  and  many  materials,  thermal  conductivity  increases  v/ith 
temperature  so  the  largest  temperature  gradients  will  be  near  the  cold  end  of  the  stack.  The  second 
reason  that  the  assumed  simple  form  for  the  temperature  distribution  is  invalid  occurs  for  the  prime 
iTtOver  above  onset  The  heat  flow  eq'iation  (20)  has  a  term  proportional  to  temperature  gradient 
and  dependent  on  the  product  of  pressure  and  paiticie  velocity.  This  term  is  a  dynamic  coefficient 
,  of  thermal  conductivity  that  depends  on  position  and  the  square  of  the  pressure  amplitude. 

For  a  given  heat  supply  the  prime  mover  can  reach  a  steady  equilibrium  state.  Heat  supplied 
at  ihe  hoc  end  is  converted  into  acousdc  energy  and  also  transported  down  the  temperature  gradient 
to  the  cold  end.  No  heat  is  transferred  from  the  stack  to  an  external  source  or  sink.  Consequently 
the  energy  flow  on  the  left  hand  side  of  Eq.  (5)  H2  is  a  constant  in  the  stack.  Merkli  and 
Thomaim,  Ref.  6,  showed  that  dH2/dz  is  proportional  to  the  amount  of  heat  transfeinng  to  the  wall 
at  location  z,  and  the  net  heat  transfering  to  the  wail  is  naturally  the  integral  over  the  same  the  wall 
length  in  question.  Swift^’^  u^  the  fact  of  H2  constant  to  determine  numerically  the  steady  state 
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equilibrium  of  the  prime  mover.  The  prime  mover  can  be  analyzed  from  a  steady  state  perspective 
as  fellows.  Given  boundary  or  iidtial  conditions  are:  heat  supplied  at  the  hot  exchanger,  the  cold 
end  temperature,  and  the  ambient  pressure  in  the  tube.  Determine:  The  acoustic  pressure, 
temperature  at  the  hot  end,  and  the  frequency  of  operation.  Solution:  Assume  an  acoustic  pressure 
at  the  rigid  end.  of  the  prime  mover.  Assume  a  hot  end  temperature  and  frequency  of  operation. 
Use  impedance  translation  to  get  the  pressure  and  impedance  at  the  hot  heat  exchanger,  stack 
interface.  Compute  H2  =  heat  supplied  at  the  hot  exchanger  plus  H2  computed  from  use  of 
equation  (5).  Since  H2  is  a  constant  in  the  stack,  the  value  computed  at  the  hot  neat  exchanger, 
stack  interface  is  the  same  throughout  the  stack.  Solve  £q.  (5)  for  the  temperature  gradient  and 
numerically  integrate  impedance,  pressure,  and  temperature  through  the  stack  to  the  cold  end. 
Here  you  must  cid  up  with  the  temperature  being  equal  to  the  cold  end  temperature:  if  not,  then 
some  of  the  assumed  quantities  are  wrong.  Assuming  that  the  loop  has  determined  the  correct 
pressure  and  hot  end  temperature,  the  specific  acoustic  impedance  looking  into  the  cold  heat 
exchanger  must  be  equal  10  the  SAI  looking  down  the  tube  towards  the  tennination  at  the  cold  end. 
Tf  not,  then  the  frequency  of  opererion  must  be  adjusted  and  the  code  should  start  the  whole  loop 
over.  That  impedance  looking  each  way  is  the  same  simply  occurs  due  to  conservation  of  pressure 
and  volume  velocity  in  the  tube.  An  alternative  approach,  or  at  least  a  mechanism  to  use  to  check 
for  energy  conservation,  which  has  some  merit  is  energy  conservation.  The  heat  entering  the 
system  must  be  equal  to  the  sum  of  heat  leaving  at  the  cold  end  plus  acoustic  energy  by  rhe  first 
law  of  thermodynamics. 

SUNDERY  REMARKS 

A  series  of  interesting  fundamental  investigations  have  been  performed  using  a  single  circular 
tube  for  which  the  thermal  boundary  layer  thickness  was  approximately  equal  to  the  tube  radius.^' 
31  Their  basic  arrangement  was  a  tube  closed  at  one  end  by  a  pressure  transducer  and  at  the  other 
by  varying  transducers  depending  on  the  particular  situation  under  study.  The  gaseous  helium 
filled  tube  was  bent  to  a  U*shape  and  the  U*portion  was  immersed  in  cold  helium  gas  at 


temperatures  in  the  range  4,2  to  45  K.“^  The  hot  end,  was  held  at  room  temperature.  They 
iiivestigated  the  stability  curve  of  the  second  tube  mode,^'  the  stability  curve  of  the  fundamental 
mode  under  a  variety  of  conditions,^®  the  universal  properties  of  a  thermoacoustic  oscillator  at  the 
intersection  of  the  stability  curve  of  the  first  and  second  modes,^^  and  the  universal  properties  of  a 
driven  thermoacoustic  oscillator. On  the  practical  side,  thermo  acoustics  refrigeration  has  been  a 
recent  topic  of  investigation. ^2-35  Traveling  wave  thermoacoustic  engines  utilizing  the  Stirling 
thermodynamic  cycle  have  also  been  of  recent  interest 

Future  work  might  include  the  following.  Numerically  simulate  the  Lagrangian  picture  of 
thermoacoustics  including  viscosity  and  temperature  gradients  and  a  variety  of  capillary  tube 
shapes.  This  would  enhance  the  intuitive  understanding  of  theimoacoustics  and  would  possibly 
lead  to  new  designs  for  thermoacoustic  elements.  Construct  and  analyze  a  prime  mover  where  the 
tubes  have  increasing  radii  from  the  cold  to  hot  end  so  that  the  optimal  thermal  boundary  layer 
thickness  (which  is  temperature  dependent)  is  achieved  at  all  positions.  Theory  indicates  variable 
tube  radii  should  lead  to  more  efficient  operation.^-3*  A  major  effort  could  be  expended  in  both 
theoretically  and  experimentally  determining  the  role  of  acoustic  streaming^  in  thermoacoustics. 
The  present  theory^^  of  theimoacoustic  streaming  is  only  a  boundary  layer  theory.  Someone 
should  build  and  investigate  more  thoroughly  the  radial-wave  thermoacoustic  engine.^  A 
comprehensive  design  program  could  and  should  be  written  for  thermoaccustic  engines,  given  the 
present  state  of  theoretical  de'.  elopment  Stability  curves  like  the  type  described  in  Refs  28-31 
could  be  determined  for  arraiigements  more  typical  of  practical  devices  such  as  that  shown  in  Fig. 
1.  We  have  recently  completed  construction  of  a  minimal  ATonset  closed  heliuin  prime  mover 
system  which  is  fully  instrumented  md  will  be  analyzed  in  the  near  future.  Much  of  the  numerical 
work  for  this  system  has  already  been  done. 

Swift  has  recently  suppFed  me  with  a  draft  of  a  paper  on  the  performarwe  of  a  large 
thermoacoustic  prime  mover.  In  it  he  establishes  a  first  crack  at  the  mature  stage  of 
theimoacoustics.  Though  we  all  know  nonlinearity  is  a  general  frontier  at  the  moment.  Swift  and 
to  some  extent  Muller^  were  first  to  establish  calculation  techniques  in  theimoacoustics.  As  Ralph 


Gcxxlman  once  said,  and  I  paraphrase,  linear  acoustics  is  women’s  work.  He  was  paraphrasing  an 
Englishman  himself.  However,  after  viewing  the  work  of  women  at  the  Univ.  of  Texas  on 
nonlinear  acoustics,  I  think  his  statement  should  have' "women’s”  replaced  with  "wimp's”.  One 
obvious  feature  of  prime  movers  is  that  the  large  acoustic  amplitudes  result  in  harmonic 
production.  Swift  treated  the  1st  harmonic  as  being  due  to  a  source  at  the  end  and  included  the 
thennoacoustic  effects  in  a  calculation.  I  see  no  reason,  other  than  wimpdom,  to  quit  at  the  first 
harmonic.  One  must  first  learn  how  to  compute  the  amount  of  harmonic  produced  for  a  given 
fundamental.  Then  with  a  stead  fast  relation  between  harmonics  and  the  fundamental,  the  overall 
system  can  be  treated  much  in  the  same  manner  as  approached  by  Swift  He  neglected  the  effects 
of  the  acoustic  wind  (or  streaming)  in  his  calculation,  but  some  accounting  of  the  energy  consumed 
by  it  should  be  possible.  Though  my  statements  here  arc  very  sketchy,  Swift's  work^^  indicates 
this  is  an  approachable  frontier  area  of  thermoacousdes.  He  finds  satisfactory  agreement  in  almost 
all  cases  when  linear  thermoacoustics  is  valid.  It  is  necessary  to  understand  the  nonlinear  problem 
if  one  wants  to  design  practical  high  power  theimoacoustic  engines. 
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FIGURE  CAPTIONS 

Fig.  1.  a)  Generic  arrangement  used  in  thermoacoustic  heat  engines,  b)  An  exposed  view  of  a 
thermoacoustic  element  consisting  of  a  parallel  combination  of  square  capillary  tubes, 
c)  A  single  arbitrary-perimeter  capillary  tube  for  use  in  a  thermoacoustic  element 


Fig.  2.  a)  Parallel  plate,  b)  circular,  c)  rectangular,  and  d)  equilateral  triangular  capillary  tube 
geometries. 

Fig.  3.  Lagrangian  view  of  a  fluid  parcel  in  a  standing  wave  near  a  boundary. 

Fig.  4.  Magnitude  and  phase  of  the  excess  temperature  between  parallel  plates. 

Fig.  5.  Real  and  imaginary  part  of  the  F(X)  different  pore  geometries. 

Fig.  6.  Demonstration  thermoacoustic  oscillator  and  analysis  impedance  tube. 

Fi^.  7.  Real  part  of  the  specific  acoustic  impedance  at  the  mouth  of  the  prime  mover. 

Fig.  8.  Imaginary  part  of  the  specific  acoustic  impedance. 

Fig.  9.  Prime  mover  sound  spectrum  for  onset  AT  =  176  K  and  a  higher  AT  =  209  K. 

Fig.  10.  Prime  mover  sound  spectrum  for  AT  =  209  K  in  the  mouth  and  5.7  cm  inside  of  the 
prime  mover. 

Fig.  11.  Spectral  peaks  as  a  function  of  distance  from  the  mouth  of  the  prime  mover. 

Fig.  12.  Time  evolution  of  the  superheated  prime  mover. 

Fig.  13.  Scale  drawing  of  the  constructed  UM  helium  filled  prime  mover. 

Fig.  14.  Stability  curves  for  the  fundamental  and  fint  harmonic  for  the  UMTAE. 

Fig.  15.  a)  Quality  factor  and  b)  resonant  frequency  for  the  UM  TAE  with  ambient  pressu'  e  173 
kPa 
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■  APPENDIX:  FORTRAN  PROGRAM  FOR  ANALYSIS  OF  THE  UM  TAE 

EXECUTION  FILE  FOR  THE  THEHMOACOUSTTC  ENGINE  CODE.  THE  CODI 
IS  ON  MY  ACCOUNT  ON  THE  IBM  3084.  CMS  OPERATING  SYSTEM 
MACHINE.  ACCESS:  LOGIN  PAARNOTT.  PASSWORD  FROGGY. 

&IF  &INDEX  NE  0  AGOTO -HEREWEGO 

&TYPE  ENTER  THE  NAME  OFTHE  PARAMS  FILE  WHICH  DESCRIBESTHETAE. 

&READARGS 
-HEREWEGO 
VMFCLIiAR 

COPYFILE  Al  PARAMS  A  TEMP  PARAMS  A  (REPLACE 
CiECGLOBALS 
HILEDEF 1 DISKT  DATA  A 
FILEDEF2  DISK  TEMP  PARAMS  A 
FILEDEF 13  DISK  PRESS  DATA  A  (RECEM  V 
FILEDEF 16  DISK  WORKFLOW  DATA  A  (RECFMV 
FILEDEF 17  DISK  HEATFLOW  DATA  A  (RECFM  V 
FILEDEF  18  DISK  IMPED3D  DATA  A  (RECFMV 
FILEDEF  19  DISK  IMPED  DATA  A  (RECFMV 
FILEDEP20  DISK  ENTHALPY  DATA  A  (RECFMV 
FILEDEF21  DISK  WORKDRIV  DATA  A  (RECFM  V 
FILEDB'27DISKQUALREC  DATAA 
FILEDEF28DISXRESFREQ  DATAA 
FILEDEF29DISKQUALFACr  DATAA 
FILEDEF30  DISK  POSVSnT  DATA  A 
FILEDB>37DISKQUALRH:  DATAEIGA 
FILEDEF  38  DISK  RESFREQ  DATAEIG  A 
FILEDEF39 DISK QUALFACr  DATAEIGA 
FILEDEF41DISKTQZ  DATAA 
FILEDEF  42  DISK  REIMFED  DATA  A  (RECFM  V 
FILEDEF  43  DISK  IMIMPED  DATA  A  (RECFM  V 
FILEDEF  44  DISK  STABCURV  DATA  A  (RECFM  V 

•  PROGRAM  FOR  COMPUTING  PROPAGATION  CONSTANTS-. 

«LOADTAEimL(  CLEAR  NCMAP  START 

•  GENERAL  PROGRAM  FOR  COMPUTING  ZPJIEAT.WORK. 

•LOADTAQ1E2  (CLEAR  NOMAP  START 

•  GQIERAL  PROGRAM,  FINITE  DOT  IN  STACK.  FOR  CCMPUTINC  Z,PJIEAT.WORK. 

•LOADTAEKE3  (CLEAR  NOMAP  START 

•  GENERAL  PROGRAM,  FINITE  DOT  IN  STACK,  FOR  COMPUTING  iPJHEAT, WORK. 

•  GAS  ABSORPTION  IS  INCLUDED  INTHE  OPENTUBE  SECTIONS. 

•LOADTAEHE3V2  ( CLEAR  NOM.AP  START 

•  GENERAL  PROGRAM,  RUNGE  KTHTA  IN  STACK,  FOR  COMPUTING  ZFJIEATWORK. 

*LOAD  TAERUNGE  (  CLEAR  NOMAP  START 

•  GENERALPROGRAH  RUNGE  KUTTA  IN  STACK,  FOR  COMPUTING  ZFJIEAT.WORK. 

•  CAN  ALSOCOMPUIETRAVEUNG  WAVES.  SOMEWHATOFIIMIZEIX 
•LQADTAEHEUNG  ( CLEAR  NOMAP  START 

•  GENERAL  PROGRAM  FOR  COMPUTING  ZJ>,HEAT.WORK.  BASED  ON  DF/DLAMBDA. 

•LOADTAEHE4  (CLEAR  NOMAP  START 

•  PROGRAM  FOR  COMPUTING  TTIE  Q  FROM  COMPLEX  EICENFREQUENCY. 

*LOADTAEHBQ  (CLEARNOMAPSTART 

•  PROGRAM  PC»  COMPUTING  THE  Q  FROM  COMPLEX  EIGENFREQUENCY. 

•  USES  INPUT  FROMTAEHE3  AS  INFUTtO  STARTTHINCS  OFF  WITH. 

•LOAD  TAEHBQ2  ( CLEAR  NOMAP  START 

•  ALL  RUNGE-KUTTA  PROGRAMS  ABOVE  HERE  ARE  NOT  CORRECTED. 

•  PROGRAM  FOR  a»1PUnNO  THE  Q  AS  A  FUNCTION  OFTEMPERATURE  IN  AN  EASY 
•WAY.  USES  RUNGE  KUTTA  INTEGRATION. 

•  SEC(»ID  PROGRAM  IS  A  COMPLEX  NUMBER  VERSION  TO  GETTHE  COMPLEX  EIG  FRE. 

•LOAD  TAEAUTO  ( CLEAR  NOMAP  START 
•LOADTAECAUrO(  CLEAR  NOMAPSTART 

•  VERSION  OFTAECAUTO  FOR  AIR  FILLEDTUBES. 

•  LOADTAECAIR  ( CLEAR  NOMAP  START 

•  TO  PLOT  LQ  AND  RES  VERSUS  LAMBDA.  USING  TAECAIR.  FOR  A  PRIMEMOVER. 

•  LOADTAECBALT(  CLEAR  NOMAPSTART 

•TAECMONT:  PROGRAM  BASED ONTAECAUTOTOCOMPUIlTHEQ  AND  RESFREQ  FOR 
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•THEMONTEREYTUBE.  TAKES  INTO  ACCOUKT  THE  DEPENDENCE  ON  TEMPERATURE 

•  WHEN  (XIMPITONG  THE  THERMAL  CONDUCnviTY  OFTHE  STACK. 

•LOAD  TAECMONT  ( CLEAR  NOMAP  START 

•  UMTAE  FORTRAN :  PROGRAM  FOR  DESIGN  OF  THE  UMTAE.  USED  TAECAUTO  AS  THE 

•  STARTING  PROGRAM. 

•LOAD  UMTAE  ( CLEAR  NOMAP  START 

•  UMTAECH:  CHECK  OF  UMTAE  WITH  MONTEREY  DATA 
•LOAD  UMTAECH  (CLEAR  NOMAP  START 

•  PROGRAM  FROM  TAECBALT  FOR  COMPUTINr  THE  RESPONSE  OFTHE  DEMOTAE  IN  AIR 
•LOAD  TAEAIRD  (CLEAR  NOMAP  START 

•  UMrAEV2  FOIOHRAN ;  PROG  FOR  DESIGN  OFTHE  UMTAE  USED  TAECAUTO  ASTHE 
•STARTING  PROGRAM. 

LOAD  UMTAEVl  ( CLEAR  NOMAP  START 

ERASE  TEMP  PARAMS  A 

ERASE  FILE  SCRATCH  A 

FUEDEF*  CLEAR 

cSEXIT 

TAB  PARAMETER  FILE.  ACTUAL  #S  OF  THE  UMTAE  TUBE. 

TAE  PARAMETER  RLE  ACTUAL  #S  OF  THE  UMTAE  TUBE 

DEFINETHETUBE  FROM  RIGHTTO  LEFT.  AT  LEFT  ONE  USUALLY  HAS  THE  DRIVER. 
MINIMUM  FREQUENCY,  MAXIMUM  FREQUENCY,  (HZ),  AND  NUMBER  OF  FREQ.  POINTS. 
2SO.OOO  400.  300 

TERMINATION  ATTHE  RIGHT  ENDOFTHETUBE 
ONE  OF  RIGID,  FREE  OR  INFIN.  INFIN  IS  AN  INFINTTEnJBE 
RIGID 

AMBIENTPRESSUREINTHETUBE  ANDTHE  DRIVER  PRESSURE  AMPLITUDE  FOR  ALL 
FREQUENCIES.  PRESSURE  IN  PASCAL,  DRIVER  DISPLACEMENTIN  METERS. 

3.0003  l.D>< 

NUMBER  OF  SECTIONS  IN  THE  TAE  EC.  AN  OPEN  SECnON,  RGH  HEAT  EXCH, 
STAOLLEFTHEATEXCH,  AND  ANOTHER  OPEN  SECnONWOULD  BBS.  INTEGER. 

5 

definition  op  SECTION  1  •«••••••••••••••••*•••••«••••••••••*• 

SECnONTYPE,  ONBOFOPENTU,  HEXCH  OR  STACK. 

OPENTU 

ELEMENTTYPE  HAS  MEANING  ONLY  FOR  HEXCH  OR  STACK  SECTION  TYPES. 

ONE  OF  RECT,  CYL,  OR  SUT,  DEFINING  THE  TYPE  OF  PORES. 

SUT 

NUMBER  OF  LAYERS  THIS  SECnON  IS  BROKEN  UP  INTO. 

1<»  NUMLAY  <»  100  PRACTICALLY. 

1 

LENCTHOFTHESECnON. 

METERS. 

23.07D-2 

TEMPERATURE  OFTHE  RGH  END  OFTHE  SECTION.  FOR  AN  ISCTHERMALSECnON 
SUCH  AS  OPEN  TUBE  OR  HEAT  EXCHANGERS,  USETRCH«n.EFT.  KELVIN. 

TEMPE^TUREATTHELEFTENDOFTkESECnON. 

SEENOTEABOVE.  KELVIN. 

293. 

RATIO  OF  2  FORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYURADIUS,  SLTT- 
WIDTH,  R£CT>2SW  A/(I>A)  A>  1>SIDES  ASPECT  RATIO,  SW-SHORTEST  SEMIWIDTH. 
4SS1D-2 

ASPECTRATIOOFTHEPORE.VALIDPORRECTANCUIARPORESONLY. 

NECESSARY  ASACENERALRULE. 

1.0DO 

POROSTTY  OF  THE  SECTION.  FOR  OPEN  TUBE,  USE  POROSITY  « I. 

FOR  OTHER  TYPES  OF  SECTIONS,  POROSTTY  <■!. 

IJX) 

END  OP  SECTION  L 

DEFINTnON  OP  SECTION  2 
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SECTION  TYPE,  ONE  OF  OPENTU.  HEXCH,  OR  STACK. 

HEXCH 

ELEMENTTYPE.  HAS  MEANING  ONLY  FOR  HEXCH  OR  STACK  SECTION  TYPES. 

ONE  OF  RECT,  CYT,  OR  SUT,  DERNING  THE  TYPE  OF  PORES. 

SUT 

NUMBER  OF  LA YERSTHIS  SECTION  IS  BROKEN  UP  INTO. 
lc=  NUMLAY  <=  100  PRACTICALLY. 

1 

LENGTHOFTHESECnON. 

METERS. 

M14D-2 

TEMPERATURE  OFTHE  RGH  END  OFTHE  SECTION.  FOR  AN  ISOTHERMAL  SECTTION 
SUCH  AS  OPEN  TUBE  OR  HEAT  EXCHANGERS.  USETRGH = TLEFT.  KELVIN. 

293. 

TEMPERATUREATTHELEFTENDOFTHESECnON. 

SEE  NOTE  ABOVE.  KELVIN. 

293. 

RATIO  OF  2  PORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYURADIUS,  SUT> 
WIDTH.  RECT>2SW  A/(1>A)  A^lsSIDES  ASPECT  RATIO.  SWaSHORTEST  SEMIWIDTH. 
1.118D-3 

ASPECT  RATIO  OFTHE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY. 

NECESSARY  AS  A  GENERAL  RULE.  ASPECT  RATIO  «>  1  always. 

1.0D0 

POROSmr  OFTHE  SECTION.  FOR  OPEN  TUBE,  USE  POROSITY  » I. 

FOR  OTHER  TYPES  OF  SECTIONS.  POROSITY  <xl. 

.(S4D0 

END  OF  SECTION  Z 

•••••••«,  definttion  of  section  3 

SECTION  TYPE,  ONE  OF  OPENTU.  HEXCH.  OR  STACK. 

STACK 

ELEMENTTYPE.  HAS  MEANING  ONLY  FOR  HEXCH  OR  STACK  SECTION  TYPES. 

ONE  OP  RECT,  CYU  OR  SLTT.  DEFININa  THE  TYPE  OF  PORES. 

RECT 

NUMBER  OF  LAYERS  THIS  SECTION  IS  BROKEN  UP  INTO, 
loi  NUMLAY  <a  100  PRACTICALLY.  (WAS  20  AT  ONE  TIME) 

2S  , 

LENCTHOFTHESECnON. 

METERS. 

5.0SIV-2 

TEMPERATUREOFTHERCHENDOFTHESECnON.  FOR  AN  ISOTHERMAL  SECTION 
SUCHAS0PENTUBE0RHEATEXCHANGERS,USETRGH«TL£FT.  KELVIN. 

293. 

TEMPERATUREATTHELEFTENDOFTHESECnON. 

SEENOIEABOVE  KELVIN. 

293. 

RATIO  (V  2  PORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYURADIUS.  SUT« 
WIDTH.  RECT«2SW  A/(UA)  A>1«SIDES  ASPECT  RATIO.  SW«SHORTEST  SEMIWIDIH. 
.T7D-3 

ASPECT  RATIO  OFTHE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY. 
NECESSARY  AS  A  GENERAL  RULE.  ASPECT  RATIO  » 1  ALWAYS. 

LODO 

POROSnr  OFTHE  SECTION.  FOR  OPtN  TUBE,  USE  POROSITY  ■  I. 

FOR  OTHER  TYPES  OF  SECTIONS,  POROSITY  <■!. 

.6900 

END  OP  SECTION  3. 

••••••••  definition  of  SECTION  4 

SECnONTYPE,  ONE  OPOPENTU.  HEXCH,  OR  STACK. 

HEXCH 
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ELEMENTTYPE.  HASMEANINGONLYFORHEXCHORSTACKSECTIONTYPES. 

ONE  OF  REGT,  CYL,  OR  SUT,  DEFINING  THE  TYPE  OF  PORES. 

SUT 

NUMBER  OF  LAYERS  THIS  SECTION  IS  BROKEN  UP  INTO. 

!<=  NUMLAY  <=  100  PRACTICALLY. 

1 

LENGTH  OFTHESECTION. 

METERS. 

1.638D-2 

TEMPERATU  EOFTHERGK  END  OFTHESECTION.  FOR  AN  ISOTHERMAL  SECTION 
SUCHASOPENTUBEORHEATEXCHANGERS,USETRGH=TLEFT.  KELVIN. 

293. 

TEMPERATUREATTHELEFTENDOFTHESECnON. 

SEE  NOTE  ABOVE.  KELVIN. 

293. 

RATIO  OF  2  PORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYURADIUS,  SUT= 
WIDTH.  RECT=2SW  A/(l+A)  A>1=SIDES  ASPECT  RATIO,  SW=SHORTESTSEMrWlDTH. 
1.118D-3 

ASPECT  RATIO  OFTHE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY. 
NECESSARYAS  AGENERALRULE.  ASPECTRATIO=>  I  always. 

1.0D0 

POROSTTY  OFTHE  SECHON.  FOR  OPEN  TUBE.  USE  POROSITY  =  1. 

FOR  OTHER  TYPES  OP  SECTIONS,  POROSITY  0.1. 

:640D0 

END  OF  SECnON  4. 

.,•••••«  DEFINITION  OF  SECTION  5 
SECTION  TYPE,  ONE  OF  OPENTU,  HEXCH,  OR  STACK. 

OPENTU 

ELEMENTTYPE.  HAS  MEANING  ONLY  FOR  HEXCH  OR  STACK  SECTION  TYPES. 

ONE  OF  RECr,  CYL,  OR  SITT,  DEFINING  THE  TYPE  OP  PORES. 

SUT 

NUMBER  OF  LAYERS  THIS  SECTION  IS  BROKEN  UP  INTO, 
loi  NUMLAY  o  100  PRACnCALLY. 

1 

LENGTHOPTKESECnON. 

METERS. 

1.29 

TEMPERATUREOFTHERGH  END  OFTHE  SECTION.  FOR  AN  ISOTHERMAL  SECTION 
SUCH ASOPENTUBEORHEATEXCHANGERS, USETRGHaTLEFT.  KELVIN. 

293. 

TEMPERATUREATTHELEFTENDOPTHESECnON. 

SEE  NOTE  ABOVE  KELVIN. 

293. 

RATIO  OF  2  PORE  AREA  TO  PORE  PERIMETER  (M).  FOR:  CYURADIUS,  SUT- 
WIDTH,  RECr>aSWA/(UA)  A>t«SlDES  ASPECTRATIO,  SW-SHORTESTSEMIWIDTH 
4:Z61D-2 

ASPECT  RATIO  OFTHE  PORE.  VALID  FOR  RECTANGULAR  PORES  ONLY. 
NECESSARY  AS  A  GENERAL  RULE.  ASPECT  RATIO  1  ALWAYS. 

1.0D0 

POROSTTY  OFTHE  SECTION.  FOR  OPENTUBE,  USE  POROSITY  ■  1. 

FOR  OTHER  TYPES  OP  SECTIONS,  POROSITY  ol. 

.995D0 

END  OF  SECTION  5. 
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The  main  program  is  the  only  part  of  this  program  that  is  specific  to  the  UM  TAB. 

••SUBROUTINE  SETVAL  •  GETTHE  INPUT  PARAMETERS  FROM  AN  EXTERNAL  FILE** 
SUBROUTINE  SETVAL 
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•  variables  WHICH  DEFINE THETAE 
CHARACTER*70SECTYP(100).ETYPE(100).TERM!N 
INTEGER  NUMl-\Y(100XNUMSEC.NUMFRE 
RiiAL'8  DELEM(100),TRGH(  100).TLEFT(  100)JIAT10(100) 

•^PEC  r(  100).POROS(  100), 

•  THCOND(100)4IECAP(100).PAMB.FREMIN.FREMAXDDR1VE 
INTEGERJ 

CHARACTER  DUMMY 
COMMON /VARSl/ ScCTYP,ETYPE,TERMIN 
COMMON /VARS2/ NUMLAYJ^MSECNUMFRE 
COMMON  /VARSV  DELEM.TRGH.TLEFT.RATIO, 

•  MPECTJ>0R0S,THC0NDJIECAP.PAMB,FREMIN,FREMAX45DR1VE 

1  FORMAT  (Al) 

2  FORMAT  (A70) 

REWTND2 

READ  (Zl)  DUMMY 
READ  (2,1)  DUMMY 
READ(Z1)DUMMY 

read  (2,*)  FREMIN.FREMAXJ^in^FRE 
READ  (2.1)  DUMMY 
READ  (2.1)  DUMMY 
R£AD(2Z)'nERMIN 
CALLNOPADCTERMIN) 

READ  (2,1)  DUMMY 
READaDDUMMY 

READ  a*)  PAMBJDDRIVE 
READaDDUMMY 
READaDDUMMY 
READa*)NDMSEC 

DO  lOJ-NUMSEQl.-l 
READ  (24)  DUMMY 
READaDDUMMY 

READ(2J)SBCTYP(D 
CALL  NOPADCSSCTYPO)) 

READ  (24)  DUMMY 
READaDDUMMY 
READ(25)ErYPE(J) 

CAU.NOPAD(ETYPE(D) 

READaDDUMMY 

READaDDUMMY 

READa*)NUMLAY(D 

READaDDUMMY 

READaDDUMMY 

READa*)DELEM(D 

READaDDUMMY 

READaDDUMMY 

READa*)TRCH(D 

READaDDUMMY 

READaDDUMMY 

READa*)TLEnXJ) 

READ  (24)  DUMMY 
READaDDUMMY 
READa*)RAT10(D 
READaDDUMMY 
READaDDUMMY 
READa*)ASPECr(D 
READaDDUMMY 
READaDDUMMY 
READa*)POROS(D 
READaDDUMMY 
THCOND(X)-0.D0 
10  HECAP(J)^.0D0 
REW1ND2 
RETURN 
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END 


•  SUBROUTINE  NOPAD  •  GETS  RID  OF  BLANKS  IN  NAMES 


SUBROUHNE  NOPAD(WAME) 
CHARACrER*70  OLDNEW.NAME 
CHARACTER*!  0(70)N(70) 
INTEGER  U 

equivalence  (OLDAD) 
equivalence  (NEW.N(1)) 
old-name 

NEW- 

• 

J-0 

do  101-1.70 

NOV 

IF  (0(1)  NE.'-)  THEN 

J.J>t 

N(J)-0(I) 

END  IF 

10  CONTINUE 
NAME-NEW 
RETURN 
END 


•  SUBROUTINEWFLOW*OOMPUTCTHE  WORK  FLOW  ATZ 


SUBROUTINE  WFIjOW(P1ZW2) 

REAL’S  W2 
COMPLEX*16  P12 

W2-CDABS(P1)*»2  •  DIMAG((aD0.1iM)/Z>  /  ZOO 
RETURN 
END 


•  SUBROUnNEQFLOW'COMPUTETHEHEATFlOWATZ.  ••••••••• 

•  IHAVEASSUMEDTHECOEFFIOENrOFTHERMALEXPANSIONBErA-l/rEMP.** 

SUBROUTINE  QFLOW(?OROSJ>1ZFLAMJ1AMT.W.DENS,TQZJCCASKSOUD.Q2) 

REAL’S  POROSDENS,TQZJCGASJCSOUD,Q2.NPR.CP.GAMMA 

COMPl£X’I6  PIZFLAMJFIAMT.W 

COMMON  /PHYCON/  GAMMAi^PR.CP 

Q2-POROS*CDABS(P1)’CDABS(P1)/2.DO 

Qi-Q2’DIMAC((OJ».IJ)0)’(OCONJC(FLAMTVFLAM-1.DO)/ 

’  (P<»OS’Z’(IJDO-NPR))- 
’  TOZ  •  DENS  ’  CP’(FLAM’NPR  ♦DCONJC(F'-AMrr) )  / 

,  •  (POROS”2*W’(CDABS(FLAM•Z))”2’(liX>■NPR”2)))- 
’  TnZ’(TOROS*KCAS  -  (1  JX>>POROS)’KSOUD) 

RETURN 


•  SUBROUnNEDERlVS’COMPUrBTHBDERlVATlVESDZUZANDDP/C;ZPOR 
’  THE  RUNG&KUTTA  WORK. 


SUBROUTINE  DERIVS(ZErA.ZINT.ALPRIMJ»ZDPDZJDZDZ) 
COMPLEX’ 16  2ETAjaNT.ALPRIMJ»ZDPDZJDZDZiFACFAa 
I.(0.D0.1J30) 

FAC-Z/ZINT 

PAC-(1J»-FAC*PAC) 

FAa-I’ZETA’ZTNT 
DZDZ  -  PAC2  ’  FAC  ♦  IDO  ’  ALPRIM  ’  Z 
’  FAC-l’ZETA’(21NT-Z’ZOTfr) 

’  DZDZ-FAC>ZD0’ALPR1M’Z 
DPtXS-FAa’P/Z 
RETURN 
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END 


•  SUB  CTAKPM  •  GET  THE  MANY  PARAMETERS  WHICH  ARE  TEMPERATURE  DEPENDENT 

•  IN  THE  STAK.  USED  PDR  RUNCE  KUTTA  INTEGRATION. 

•  THE  STACK  IS  ASSUMED  TO  HA  VE  POROUS  WALLS 


SUBROUTINE  STAKPM(ErYPE.W,POROUSPAMB.T.TaZJ3ENSRATiaASPECT. 
•  FLAM.FLAMrZErA/.LPRIM2IND 

CHARACrER*70  ETYPE 

REAL'S  P0R0USPAMB.T,TCZJ)ENSRATI0,ASPECT 
COMPLEXM6  aAMJlAMr2ErA^RlMJINr.W.LAMBDAXAMBDT 

•  POROUS  WALL  VARIABLES 

REAL'S  PORTOrjORWAUDWALUWFACT 
a»<PLEX*16XI 

'  SUPPORTING  ROLE  VARIABLES 

REAL'S  SSPEED.VISCCPJ4PR,GAMMA 
COMMON  /PHYCON/  GAMMA.NPR,CP 
'  SETTHE  POROUS  WALL  CONSTANTS  TOR  THE  aX)  CELL  CERAMIC. 

'TOTURN  OFFTHE  POROUS  WALL  CALCULATION  JUST  LET  BELOW. 

FORWAL- 04900 
DWALL- 100  JV6 

POtTOr  -  POROUS'd  ♦  ZDO'PORWAL'DWALLfRATlO) 

WFACT  -  GAJUMA'tPORTOT-POROUSVaDO'POROUS) 

•  begin 

CALL  VDCHE(T^AMB.VISCJ}ENSSSPEEaKCAS) 

CALL  GETIAM(DENSVISC.WJUT10LAMBDA) 

LAMBOT  .  DSQRTINPR)  •  LAMBDA 
IF  (ErYPEL£Q.-RECn  THEN 
CALL  FR£CT(ASPECTLAMBOAJFlAM) 

CALL  FRECT(ASPECTLAMBDTJLAMD 
ELSE  IF  (EIYPEEQ-XLYL*)  THEN 
CALLPCYLOAMBDaJFLAM) 

callpcylclambdtjlamd 

ELaiF(ETYPEJQ.'SLrr)THEN 

CALLFSLrrOAMBOAJLAM) 

CALL  FSLrrOAMBOr JLAMT) 

nat 

STOP 

SJDIF 

CALLWNHEXOFLAMTJFLAM.WSSPEEDZETA) 

ZDVr  -  DENS'W  /  (POROUS  •  FLAM  •  ZETA) 

GAMMA  •  (OAMMA-iiXI)'FLAMT 
XI-1J30  +  WFACT/XI 
ZETA-ZETA'Xl 
ZlNT-ZINr/XI 

ALPROI .  TOZ  •  (FLAMT  /  FLAM  •  1J30)  /  (100  •  T  •  (1.00  •  NPR)) 

REnXRN 

END  . 


•  SUBROUriNEZFTRAN'DOTHEIMFEDANCETRANSLATlONTHEOREM**** 

•  •  DO  ALSOTHE  PRESSURE  TRANSLATION  THEOREM  *••• 

•  THIS  VERSION  IS  FOR  THE  HEATEXCHANCERSANDOPENTUBESECnONS 

SUBROUTINE  ZFTRAN(ZINT.SN.CSZZMDJ>1  J>1MD) 

COMPLEXM6Z2MDaNrSN.CSCr.PI.PIMDJAC 

CT-CS/SN 

FAC-22NT/Z 

ZMD  -  ZINT  •  ( CT  -  (O.DO  liW'PAC )  / 

•(FAC'Cr -(0.00.1.00)) 

PIMD  -  PI  •  ( CS  -  (0.D0.1D0)  •  FAC  •  SN ) 

RETURN 

END 


•  SUBROUTINE  ZRIGID**  IMPEDANCE  OF  A  RIGIDTERMINATION. 


SUBROUTINE  ZRlGID(DENS,SSPEEaVTSC.W2) 

REAL'S  DENS^PEED.V1SC.NPR.GAMMACP 

C0MPLEX*16Z.W,FAC 

COMMON /PHYCON/ GAMMA.NPR,CP 

FAC  -  CDSQRT(  DENS»SSPEED**2  /  (W'VISC) ) 

Z  -  (1.DOHI»)*DENS*SSPEED*FAC*DSQRT(NPRV(DSQRT(1DO)‘(GAMMA  - 
•  1.0D0)) 

RETURN 

END 


•  SUBROUTINE  WNTUBE  «  WAVENUMBERS  FOR  WAVES  IN  THE  OPEN  TUBE  PARTS. 


SUBROUTINE  WNTUBE(LAMBDA ,  W-.  SSPEED ,  K) 

REAL'S  SSPEED.GAMMA.NPROJVMBDAJ'ACI.CP 
COMPLEX' 16  K.W 

COMMON  /PH  YCON/  GAMMA J^CP 

FACl  -  ( 1.D0  +  (GAMMA  -  l.DO)  /  DSQRT(NPR) )  /  DSQRT(1D0) 

K  -  W/SSPEED  •  ( l.DO  ♦  (loX) .  1  J»)  '  FACl  /  LAf^BDA  ) 

RETURN 

END 


•  SUBROUTINE  FTUBE  •  COMIUTES  RLAMBDA)  FOR  THE  RESONANTTUBE 


SUBROUTINE  FTUBE(LAMBOAJFLAM) 

COMPLEX*16  LAMBDAJIAM 

FLAM  -  1.00 .  (LDQ.1JX)  *  OSQRT(ZCDO) /  LAMBDA 

RETURN 

END 


'  SUBROUTINE  WNHEX"  WAVENUMBERS  PC*  WAVES  IN  THE  HEAT  EXCHANCtRS. 


SUBROUTINE  WNHEXIFIAMT ,  FLAM .  W .  SSPEED .  K) 

REALt  SSPEEDOAMMAJ<PR.C7 
C0MPLEX*16  FUMTJIAMXW 
COMMON /PHYCON/ GAMMA.NPiLCP 

K  •  W/SSPEED  •  CDSQRTT  (GAMMA  •  (GAMMA  >  UXQIlAMr)  /  flam  ) 

RETURN 

END 


'  SUBROUTINE  VIXME  ••  VBCOSmr,  DENSITY,  AND  SOUND  SPEED  OF  HELIUM 
•  AS  A  FUNCTION  OFTEMPERATURE  AND  AMBIENT  PRESSURE.  ALSO THETHERMAL 
•coNDOcnvmr, 


SUBROUTINE  VDCHECTABS .  PAMB .  vise .  DENS .  SSPEEaXGAS) 

REAL'S  TABSJ>AMB.VISC.DENS,SSPEED.GAMMA,NPR,CPJCCAS 
COf.aWN /PHYCON/ GAMMAJ<PR,CP 
raNS  -  PAMB  •  4.0D-3  /  (  TABS  •  SS143D0 ) 

•  MY  EXPRESSION  FOR  VISOOSrTY. 

Vise  -  ISS7D-3  •  (TABS  /  r73.15D0)"0A567D0 

•  MY  EXPMSSION  FOR  THERMAL  CONDUCnVTTY:  NPR-CONSTANT. 

KOAS-VISC'CP/NPR 

•  swnrs  EXPRESSION  FOR  VISCOSTTY, 

•  V1SC-5.131D.7'TABS".644ID0 

•  SWIFTS  EXPRESSION  FOR  THERMAL  OONDUCTTVITY.  NPR  NOT  CONSTANT. 

•  KGAS -0.0044  •TABS".644I  DO 

•  NPR-VISC'CP/KCAS 

SSPEED  -  97HD0  •  DSQRT(TABS  /  mi5D0) 

RETURN 

END 


•  SUBROUTINE  FSLrr"' COMPUTES  P(LAMBDA)  FOR  PARALLEL  SLITS. 


SUBROUTINE  FSUTdAMBDA ,  FLAM) 


COMPLEX‘16  LAMBDAj^AM^RMluUGUHCTANHw^ 
SQRMI  -  (l.ODO .  -l.ODO)  /  DSQRT(  2.0DO ) 

AR-SQRMl*  LAMBDA 
ARGUM-CDEXP(-AR) 

AR-AR/2JX) 

CTANH  -  ( 1.DO  -  ARGUM )  /  ( l.DO  ♦  ARGUM ) 

FLAM  -  tX)D0  •  CTANH  /  AR 
RETURN  - 

END 


•  SUBROUTINE  PCYL  ••••  COMPUTES  F(LAMBDA)  FOR  CYUNDRICAL  PORES. 


St^ROUTINE  PCYIOAMBDA ,  FLAM) 

COMPLEXM6  FLAM,SQRI.CBSa)JOJl^GUM,LAMBDA 

INIECERN 

N-2 

SQRI  -  (l.ODO .  l.ODO)  <  DSQRT(  lODO ) 

ARGUM  -  SQRI  •  LAMBDA 
CALL  DCBJNS(  ARGUM .  N ,  CBS ) 

JO^CBSd) 

J1-CBS(2) 

FLAM  -  1.000 -( LODO  •  J1 )/(  ARGUM  •  JO ) 

RETURN 

END 


•  SUBROimNE  FRECT  ••••  COMPUTES  FOAMBDA)  FOR  RECTANGULAR  PORES. 

•  OPTIMIZED  BY  RANDYZAGAR.  21  NOVEMBER  1991. 


SUBROUTINE  FRECT(ASPECT.  LAMBDA.  FLAM) 

REAL**  PISQ,  ASPECT.  ASPSQ 
REAL**  XN.XM.FAC.TMN.TNHTMM 
REAL**  TrERM.SSUM 
INTEGER  M,N,SUMMAX 
COMPLEX*!*  SUK  FLAM.  YMN,  YNM,  YMM 
COMPLEX*!*  FAa,  LAMBDA.  TERM 
DATA  PISQ/9J*9*04404DQr 
DATA  FAO*J7Q22S64D>l/ 

C 

FACl  -  PISQ/ (LAMBDA  *(1J)0  4.  ASPeCD)**2 
FAaR-DREAL(FACl) 

FAai-DIMAG(FAa) 

C 

ASPSQ  •  ASPECT  *  ASPECT 
SUMMAX-S; 

SUM  .(OJXLOJIO) 

C 

OOSOM-1.SUMMAX.2 
XM-DROATTM)  , 

XM-XM*XM 

C 

TMM  -  XM  *  (ASPSQ  4. 1 JXJ) 

YMM  •  DCMPLXd  JX  r  FACII  *  TMM.  FAOR  •  TMM) 
SUM  -  SUM  4. 1 JX  /  (YMM  *  XM  *  XM) 

C 

DO  40  N  -  M42.  SUMMAX,  2 
XN-DFLOATTN) 

XN-XN*XN 

C 

TMN>(ASPSQ*XM4XN) 

YMN  -  DCMPLX(1.0D0  -  FACII  *  TMN.  FACIR  * TMN) 
C 

TNM-(ASPS0*XN4XM) 

YNM  -  DCMKXd.ODO  -  FAOI  *  TNM.  FACIR  ‘TNM) 
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mM  -  (l.i»YMN  ♦  1,D0/YT;M)  /  (XM  •  XN) 

TTESM  -TERM  •  DCONJG(TERM) 

SSUM  -  SJ-'M  •  DCOMJCkSUM) 

C 

*  GUARANTEED  ACCUR  AO'  FOR  THE  SUM  TO  EXPONENT/2.  EG.  SUM  -  SUM ERROR 

•  WHERE  1  J\EXP0NENT/2  -  E.RROR/SUM.  CXDMPARES  ON  THIS  LEVEL  OF  ACCUR  \CY 
•WTTHI'RECTO. 

IF  (TTcRM  .LT.  (SSUM*  I  .OD- 10))  GOTO  30 
C 

SUM -SL’M+ TERM 
40  CONTINUE 
30  CONITNUE 

FLAM -SUM ‘FAC 

RETURN 

END 


•  SUBROUTINE  FRECTO  ••••  COMPUTES  F(LAMBDA)  FOR  RECTANGULAR  PORES. 

•  OLD  SLOWER  VERSION. 


SUBROUTINE  FRECTCXASPECT .  LAMBDA ,  FLAM) 
REAL'S  PLASPECTwASPSQ 
REAL'S  F;CN,XMJAC 
INTEGER  M/«;5UMMAX 
CC»<PLEX'16  SUMFLAKYMNXAMBDA.FAC1 
*  SUMMAX  MUST  BE  AN  ODD  NUMBER.'!!!!*!!!!!!!!!!!*. 

PI  -  4.0DO  *  DATANd.ODO) 

FACl  -  PI  *  PI /( LAMBDA  *(  1.0DO  +  ASPECT )  )"2 

FAC-64.0D0/PI"4 

ASPSQ  -  ASPECT  *  ASPECT 

SUMMAX-51 

SUM -(0.0D0, 0.000) 

OO30M-SUMMAX.L-2 

XM-DFLOAIXM) 

XM-XM*XM 

DO40N.SUMMAX.1.-2 

XN-DFLOAT(N) 

XN-XN*XN 

YMN  -  1  J»  ♦  (OJXM  DO)  *  FACI  *  (  ASPSQ  *  XM  +  XN  ) 
40  SUM-SUM  +  1.0D0/(  XM'XN*  YMN) 

30  CONTINUE 
FLAM -SUM 'FAC 
RETURN 
END 


*  SUBROOnNEQUAFACCbMPUTESTHEQUALITYFACrORANDRESONANTFREQU. 
SUBROUTINE  QUAFAQAMP  JR£Q,QJL£SFRE) 

REAL'S  QJIESFRE.AMP(2000)  JREQ(2000)>IAXAMP4TIEHAFAMPHAF 
REAL'S  AMP2(2000).FREQ2(2000)^C(3)3L(3),BU(3),X5CAL£(3) 

REAL'S  X0UESS(3%FVALUE4’SCALE(3)JU>ARAM(7) 

INTEGER  JJRESJHALFJOJMDATJCOUNTJ^  JETART4BTYPE4PARAM(7) 
COMMON /QCALQ  AMP2411EQ2NUMOAT 
EXTERNAL  FUNCn 
N-3 

IPARAM(1)-0 

IBTYPE-0 

ISTART-0 

MAXAMP-0.D0 

JCOUNT-0 

D010)-l,2000 

IF  (AMP(J)  .CT.  MAXAMP)  THEN 
MAXAMP-AMP(J) 

RESFRE-FR£Q(0 

JRES-J 


36 


END  IF 

10  CONTINUE 

AMPHAF=MAXAMP/DSQRT(1D0) 

DO20J-U000 

IF  (AMP(J)  .GT.  AMPHAF)  THEN 
FREHAF  -  (FREQ(J)-*.FREQ(J-1))/2.D0 
Q  -  0^  •  RESFRE  /  (RESFRE  -  FREHAF) 

JHALF-J 
GOTOSO 
END  IF 

:0  CONTINUE 

30  DO  40  J-JHALF-1 JRES  +  (JRES-JHALF)  +  1 
JCOUNT-JCOUNT+ I 
AMP2(JCOUNT) .  AMP(J) 

40  FREQ2(JCOUNT)  -  FREQ(J) 

NUNCDAT-JCOUNT 

XGUESS(1)-MAXAMP 

XGUESSQ)  >  RESFRE 

XGUESS(3).Q 

DOSOJ.U 

XSCALE(J).1I» 

FSCALE(J)-1.D0 
BL(J)  -  XGUESS(J)  •  ^ 

50  BU(J)-XGUESS(J)*ZDO 

CALL  DBCONF(FUN(m.  N .  XGUESS ,  IBTYPE .  BL .  BU . 

*  XSCALE.FSCALE.IPARAM.RPARAM.XC.FVALUE) 
MAXAMP-XC(1) 

RESFRE -XCC) 

Q-XC(3) 

RETURN 

END 


SUBROUTINE  FUNCTl  FOR  IMSL  OPTIMIZATION 


SUBROUTINE  FUNCTKN .  XC ,  RMSERR) 

REAL*S  AMP(2000)41(EQ(2000)JCC(3XRMSERR.MAXSQJD.Q,F 
INTEGER  NUMDATJ 
COMMON  /QCALO  AMP JTIEQNUMDAT 
MAXSQ-XC(l)*XC(l) 

P0-XC(2) 

Q-XC(3) 

RMSERR-OIX) 

DO  lOJ-IJIUMDAT 
F-FR£Q(J) 

RMSERR  -  ( MAXSQ /  (liX)  ♦  (1D0*Q*(F-F0)/F0)»*2 ) 

•  -  AMP(J)*AMP(J))*»2  ♦  RMSERR 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  GETLAM 
Lambda  ••••••••••« 

SUBROUTINE  GETLAM(DENS .  vise .  W ,  R .  LAMBDA) 

REAL*8DENS,VISCJl 

COMPLEX*  16  IAMBDA,W 

LAMBDA  -  R  •  CDSQRT(  DENS  •  W/  V1SC ) 

RETURN 

END 


•  VERSION  3.0  FOR  HEUUM  BY  PAT  ARNOTT..  16  FEB  91  ••••••••••••••••*•• 

•  THIS  VERSION  USES  RUNGE  KUTTA  SOLUTION  PORTHE  INSIDE  OTTHE  STAK. 

•  TRANSLATION  THEOREMS  ARE  USED  IN  OPENTUBE  AND  HEAT  EXCHANGERS. 

•  WTHE  RADIAN  FREQUENCY  IS  ASSUMED  COMPLEX  EVERYWHERE 


SUBROUTINE  TAE(FMIN,FMAX^FREQJ=OESD 
COMPLEX*  16  W,LAMBDA,LAMBDT.F0FW.FPLUS,SQRI.KA 

•  GENERIC  VARIABLES  SUCH  AS  PHYSICAL  PROPERTIES  OF  THE  GAS. 

REAL»8SSPEED,VISC.CP.NPR.GAMMA,KGAS.KS0UD.ARES 

•  Z  DEPENDENT  ARRAYS.., 

PARAMETER  (N-iOOO) 

REAL*8  DENS(N).TAVE(N).T02(N).POROUS(N).DSUB(N)^OOR{N),Q2(N). 

•  W2(N) 

COMPLEX*16  ALPHA(N).K(N),COMDEN(N).FLAM(N),FLAMT(NI^(N),Pl(N) 

LOGICAL  INSrrAK(N) 

•  VARIABLES  WHICH  DEFINE THETAE 

CHARACrER*70SECrYP(100).ETYPE(100).TERMlN 

INTEGER  NUMLAY(lOO)jmiSECJ^MFREJ>LOTN 

REAL'S  DELEM(100).TRGH(l00).TLEFT(I00).RATlO(100),FRENEW 

•>VSPECr(100).POROS(100). 

•  THCOND(100).HECAP(100).PAMB.FREMIN,FREMAX,DDRIVE,PAMBTM 

•  DEFINE  SOME  GLOBAL  VARIABLES. 

REAL'S  PLTWOPl,FMIN.FMAX 
INTEGERNFREQ 

'  EXTRA  VARIABLES  NECESSARY  FOR  RUNGE  KUTTA  EVALUATION  OFTHE  PROBLEM. 
COMPLEX*  16  K1  JC2,K3  JC4>1 1,M2,M3.M4.DPDZ.DZDZ.TFUN 
COMPLEX*16ZErA0VLPRIMEINT.AR.SN.C3 
REAL'S  TNS.TNSMl.FOEST.LEQUrV.ROASESUM,OPL 
INTEGER  LUUPJLOWjnJMror.NS 
COMMON /PHYCON/ GAMMA^.CP 
COMMON /VARSl/ SECTYP.ErYPETERMIN 
COMMON /VARS2/ NUMLAYJ^JMSECNUMFRE 
COMMON  /  VARS3/  DELEM,TRGH,TLEFT.RATIO. 

•  aspect.poros,thcond.hecap.pamb.fremin,fi:emax.ddrive 

COMMON  /OUTPUT/  PlZQ2,WZZCOORJNSTAK 

*  ESTABLISH  SOME  OFTEN  USED  CONSTAf  <TS. 

PI-4J>. 'DATANaDO) 

TWOPI-ZDO'PI 

NPR-2i>0/3JX) 

GAMMA-5JX)/3J» 

EGAS  -  8S143D0  *  lOOO.ODO  /  4.0d 
CP  -  L5D0  •  8  JI43D0  /  4.D.3 
KSCXJD.0.16 

SQRI  -  (l.DO,l  JX))  /  DSQRT(ZDO) 

*  GET  DETAILS  OFTHETAE 

CAll.snVAL 

ARES  -  PI  *  RATIO(NUMSEC)"2 

*  ESITMATETHE  RESONANT  FREQUENCY  FROM  THE  LENGTH  AND  SOUND  SPEED  DIST. 

LSUM^JOO 

OPL-0.0 

DO  1  J-l  JflJMSEC 

LEQUIV- DELEMd)  •  POROS(J) 

•  LtiQUlV-DELEM(J) 

LSUM-LSUM>  LEQUIV 

SSPEED-DSQRT(GAMMA*RGAS»TLEFT(J))'(3.DO*TRCH(jyTLEFr(DV4.DO 
1  OPL-OPL+SSPEED'LEIQUIV 
POESTT  -  OPL  /  (2J»  •  LSUM**2) 

*  COUNTTHE  NUMBER  OF  BINS  USED...., 

NUMItrr-0 
DO  10  J.l^^UMSEC 
10  NUMrOT-NUMrOT>NUMlJVY(J) 

NUMIUr-NUMTXTT+1 

FMIN-FREMIN 

FMAX-FREMAX 

nfrEq-numfre 

RETURN 

*  ambient  PRESSURE  VARIATION  IS  HARDWIRED  IN  HERE 
ENTRYTAEl(WJOFWJPLUSJ>AMBTM) 
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PAMB  -  PAMBTM 

•  GCTTHE  SPECIFIC  ACOUSTIC  IMPEDANCE  AND  PRESSURE  AT  ALL  POINTS. 

•  START  ATTHE  RIGHT  AND  MOVETOTHE  LEFT. 

IF  CTERMIN  £Q.  •RIGID')  THEN 

CALLVDCHE(TRGH(NUMSEC)J>AMB.V1SC.DENS(NUMrCrD,SSPEED,KGAS) 

CALLZRIGlD(DENS(NtJMTt]fr).SSPEED.VlSC.W2(NUMTOr)) 

ELSE  IF  (TERMIN  .EQ.  •FREE')  THEN 

CALL  VDCHEfrRGH(NUMSEC)J>AMB.VISC.DENS(NUMTCrr).SSPEED,KGAS) 

CALL  GErU\M(DENS(NUMT0rr).VISC.W.RATlO(NUMSEC)MMBDA) 

CALL  WrmJBE(L\MBDA,W^PEEDJ<:(NUMrcrr)) 

KA  -  K(NUMTCrn  •  RATICKNUMSEC) 

Z(NUMTOrr)  -  SSPEED*DENS(NUMTOr)*KA*(KA/4.DO  -  (0.0D0.0.6DO)) 

ELSE  IF  CTERMIN  EQ. 'INFINITHEN 

CALL  VDCHECTRCH(NUMSEC)J»AMB,VISCDENS(NUMT0T).SSPEED,KGAS) 
CALLGETLAMCDENSCNUMTOD.VISC.W.RA'nOCNUMSEOXAMBDA) 

CALL  FTUBECLAMBDAJLAMCNUMrOfT)) 
CALLWNTUBE(LAMBDA,W^PEEDJC(NUMTaT)) 

ZCNUMTCTT)  -  DENSCNUMTOTT)  •  W  /  (FLAM(NUMTOT)  •  KCNUMTOT)) 

ELSE 

STOP 

ENDIF 

PKNUMTOrn-iJX) 

•  APPLYTHE  IMPEDANCE  AND  PRESSURE TRANaA'nON  THEOREMS  EVERYWHERE. 

•  WORK  FROM  RIGHT  TO  LEFT. 

DO40I.NUMSEC.1.-1 

JUP-0 

DOSOJ.U1J4UMSEC 
50  JUP.JUP*NUMLAY(J) 

JUP-NUMTOT-JUP-1 
JLOW  -  JUP  -  NUMUVYO)  I 
DSUBOUP)  .  DELEMO)  /  NUMLA Y(I) 

POROUS(IUP)  •  FOROSCI) 

•  IMPEDANCETRANSL4TE  FOR  THE  OPEN 'TUBE  OR  HEATEXCHANGER  SECTIONS. 

IF  (SBCTYPCD^OPENTU  .OR.  SECTYPCD^Q.'HEXCH)  'THEN 
CALL  VDCHE(TRCH(IXPAMB,VISCJ}ENS(JUP),SSPEEDJCGA5) 
TAVECJUFI-TRCHCI) 

TaZ(JUF)-0.0D0 
ALPHACJUP)  -  (0JX).0iX)) 

CALLGErLAM<DENS(JUP),VISC.WKA'nO0WAMBDA) 

U^OT-  DSQRT(NPR)  •LAMBDA 
IF  (SECTYP(D£Q.OPENTU) 'THEN 
CALL  FTUBE(IAMBDAFIAM(JUP)) 

CALL  FrUBE(lAMBDTJLAMT(JUP)) 
CALLWNTUBE(LAMBDA.W3SPE£DJC(JUP}} 

•  CALLPCYL(IAMBDAJLAM(JUP)) 

•  CALLPCYL(LAMBDTJIAMT(JUPJ) 

•  CALLWNHEX(FlAMT(A;P)JLAM(JUP).W^PEED.K(nJP)) 
•OTHSIWlSETHETUBESECnCffllSAHEATEXCHANGER.  FIND  ITS  GEOMETRY. 

ELSE  IF  (SECrYF(I)£Q.'HEXCH)  THEN 
IF  (ErYPE(l)iQ.'RECr)  THEN 
CALLFRECTCASPECraUAMBDAJFLAMOUP)) 
CALLFRECF(ASPECT(I)jAMBI7r41AMT(JUP)) 

ELSE  IF  (ETYPECDEQ.Xm.')  THEN 
CALL  PCYLOAMBDAFLAMCJUP)) 

CALL  PCYL(lAMBDrjFlAMT(JUP)) 

ELSE  IF  (ETYPECBEQ.'SLrr)  THEN 
CALL  PSUTTOAMBDA  JFLAM(JUP)) 

CALL  FSUr(LAMBDTJlAMT(JUP)) 

ELSE 

STOP 

•  anerrorontheinputofetypecdhasoccurcd. 

END  IF,  ■ 

CALLWNHEX(FlAMr(H;P)J*lAM(JUP).W,SSPEED4C(JUP)) 

ENDIF 


COMDEN(JUP)-DENS(JUPVFLAM(JUP) 

ZINT  -  COMDEN(JUP)  •  W  /  ( K(JUP)  •  POROUS(JUP) ) 

AR  -  KfJUP)  •  DSUB(JUP) 

SN-CDSIN(AR) 

CS-CIXXJS<AR) 

DO60J-JUPJLOW,-l 

INSrAK(J)..FALSE 

DENS(J)-DENS(JUP) 

TaVE(J>TAVE(JUP) 

TOZ(J)-TOZ(JUP) 

POROUS(J)-POROUS(JUP) 

ALPHA(J)-ALPHA(JUP) 

FL\M(J)  -  FLAM(JUP) 

FU^W-FIAMKJUP) 

K<J)-K(JUP) 

COMDEN(J)  -  COMDEN(JUP) 

DSUB(J)-DSUB(JUP) 

CALLZPTOAN(ZINT4N.CS^J+l)JZ(J).Pl(J+l).Pl(J)) 

60  CONTINUE 

ELSE  IF  (SECTYP(I)mSTACK')  THEN 

•  IMPEDANCE  TRANSLATE  FOR  THE  STACK  SECTIONS. 

TOZ(IUP)-(TRGH(D  -  TLEFT(I))  /  DELEMfl) 

NS-0 

DO70J-JUPJLOW,.l 

INSrAK<J)-.TRUE. 

NS.NS-fl 

TNS  -  TRGHO)  -  (TRCHO)  -  TLEFTO))  *  DFIjOAT(NS) 

•  /DFLOATTNUMLAYO)) 

TNSM!  -  TRGH(D  -  <TRGH(I)  -  TLEFTO))  *  DFLOATINS-l) 

•  /DFL0ATXNUMLAY(D) 

TAVE(J)  -  (TNS + TNSMI)  /  ZDO 
FOROUS(J)  -  POROUSdUP) 

0SU6(J)  -  DSUB(JUF) 

TQZ(J)-TQZ(IUF) 

•  STARTTHERUNG&ICUrTA 

CALL  SrAKPM(ErYPE(D.WJ>OROUS(JXPAMB.TNSMl.TOZ(J)J)ENS(J), 
2  RATIO(I)w^SPECT(I)41AM(I)JFLAMTXJ)2ETAw\LPRlMJI^ 
CAU.OERIVS(2^^A2IN^,ALPRIMJ>l(J•^l).Z(J••>l)J>POZJ^ZDZ)  . 

K1  — DSUB(J)*DPCZ 
Ml --DSUB(J)*DZDZ 
P1(J)  -  P1(J>1)  -f  K1/2J)0 
Z(J)  -Z(J-«-l)>Ml/2JD0 

CALL  STAKPM(ETYPE(D.W4>OROUS(J)  J»AMB.TAVE(J).T0Z(J)J3ENS(J), 

2  RATIO(I)^PE(n‘(I)JTJkM(J)JFL^MT(J)2ETAALITUMiI^ 

CALLDE!m^S(ZErA;aNr,ALPRlMj'l(JU(J)J5PDQU»3>5 

K2— DSUB(J)*DPCZ 

M2— DSUB(J)*DZDZ 

Pl(J)-Pl(J4l)-hK2/lD0 

Zq)  -Z(J4-l)4-M2aD0 

CALL  DERlVS(2E^A;ZI^^'wUJ1UMJ»l(J)AiU)PIXU3^ 

K3— DSUB(J)*DPDZ 
M3— DSUB(J)*DZLZ 
Pl())-Pl(I-*-l)>K3 
Z(J)  .Za-^1)-^M3 

CALL  srAKPM(ErYPE(D.WJ^ROUS(J)4»AMB,TNS.TOZ(J)JDENS(J). 

2  RATIO(I)^PECr(I).FLAM(J)JLAMT(J)JETAALPRIMjaNT) 

CALL  OERrVS(ZErA;DNTALPRIMJ>l(J).Z(J)J}PDZJDZDZ) 

IC4— DSUB(J)*DPOZ 
M4.-DSUB(J)*DZEZ 

Pl(J)  -  P1(J*1)  >  (K1+2J>)*K2+ZD0*K3+K4)/6.D0 
Z(J)  -  Z(J+1)  >  (MU2JX)*M2+2JX)*M3+M4V6.D0 
CX}MDENW-DENS(iVFlAM(J) 

ALPHA(I)-ALPRIM 
70  CONTINUE 
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ELSE 

•  AN  ERROR  ON  INPUT  OF  SECTYP  HAS  OCCURED. 

STOP 

ENDIF 

40  CONTINUE 

•  THE  IMPEDANCE  IS  NOW  KNOWN  AT  ALL  SPOTS  IN  THE  TAE 

•  GETTHE  ACOUSTIC  PRESSURE.  ZODOR,  WORK  AND  HEAT  FLUXES. 

•  ZCOOR(l)  -  O.ODO 

IF  (TERMIN.EQ.lNnNTrHEN 
TFUN  -  IDO/PKD 
ELSE 

TFUN  -  (0.D0,-1.D0)  •  W  •  2(11  •  DDR^^/  Pl(l) 

END  IF  !  TERMIN.EQ.TNFIN’  CONOmONAL. 

P!(1)-TFUN*P1(1) 

•  CALLWFLOW(Pl(l)J(l).W2(l)) 

DO80J-lJ^UMTOr-l 

PI(J+1)-P1(J+1)*TFUN 

•  ZCCOR(J+l) .  ZCOOR(J)  +  DSUB(J) 

•  CALLWFLOW(Pl(J+I)J<J+l).W2(J+l)) 

•  CALLQFLOW(POROUS(J).Pl(J).Z(J).FLAM(J).FLAMr(J).W. 

•  •  DENS(J).T(JZ(J)JCCASJCSOLID.Q2(J)) 

80  CONTINUE 

•  J-NUMTQT-l 

•  CALL  QFLOW(POROUS(J).Pl(J-t-l)J<J>l)JLAM{J).FLAMT(J).W, 

•  •  DENS(J).TOZ(J)JCGASJCSOUD.<32(NUMror)) 

•  D08lJ-lJ<UMTOr 

•  Q2(J)-Q2(;)*ARES 
•81  W2(j)-w:(j)«  ares 

•  THE  SIGN  OF2(l)  WAS  CHANCED  ON  8-12-91.  THE  IMPEDANCE  LOOKING  TO 
•THE  RIGHT  HASTO  BE  BQUALTO  MINUS  THE  IMPEDANCE  LOOKING  TO  THE  LEFT 

•  AND  THE  MENUS  SIGN  COMES  FROM  THE  DIRECnON  OF  PARTIdE  VELOCITY. 

POFW  .  SQRI  •  CDSQRT(DENS(l)*rS  •  SSPEED**4  •  NPR  /  (W»V1SQ  y 
•  (GAMMA- liX))  ♦  Z(l) 

FPtUS  «  FOFW  -  2D0  *  Z(l) 

RETURN 

QID 


•  SUBROUTINE  CHANCE  USED  TO  INSEXr  A  NUMBER  INTO  A  LINE  OF  A 
•SEQUENTIAL  FILE 


SUBROUTINE  CHANGE(FNUMANUM.TO) 
INTEGER  FNUMXNUMJ 
CHARACTER*80LINE 
1  PORMAT(A80) 

REAL*8TO 

OPENOJUE-SCRATCH) 

REWINDO) 

REWIND(FNUM) 

DO10J-1.S000 

READ(FNUM.l  JND-20)  LINE 

1F(INELNUM)THEN 

WRrTE(3,l)LINB 

USE 

WRrTE(3,*)TO 

ENDIF 

10  CONTINUE 
20  ENDFILE2 
REWINDO) 

REW1ND(FNUM) 

D030J-14000 
READ(3.1.£ND-40)LINB 
30  WRrTE(FNUM,l)  LINE 
40  ENDFILEFNUM 
REWINDO) 
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REW1ND(FNIIM) 

RETURN 

END 


•  PROGRAM :  EVALUATE  TAE  FOR  A  RANGE  OF  PARAMETERS. 

•  VERSION  EXPUCITELY  FOR  THE  UMTAE 

•  PAT  ARNOTT,  22  MARCH  1991,  MOD  1 M5-91. 

•  DETERMINES  THE  STABrLTTY  CURVE  FOR  THE  FIRSTTWO  MODES  AS  A  FUNCTION 

•  OF  THE  AMBIENT  PRESSURE  INCLUDES  RNTTE  WALL  POROSITY  IN  THE  STACK. 

PROGRAM  UMTAE2 

•  VARIABLES  USED  TO  GETTHE  Q. 

REAL*8  AMPLrT(2000)  JREQU(2000).QUAURESFRE4IESOLD.QOLD 
REAL’S  QOLD2,TOLD2^LD2.TOLD 

•  VARIABLES  RETURNED  FROM  TAE 

REAL’S  FREMIN,FREMAXPOEST 
INTEGER  NUMFRE 

•LOCAL  variables  TO  THE  MAINUNE....: _ 

REAL’S  PAMB,PMIN,PMAX.TMIN 

REAL’S  TRGH.TONSET^TWOP1,TEST.TESTF.TNOLD,TNOLD2 
INTEGER  TTRGHJTXnN  JFREQjmMESJIPAMBSJ^RESSJ^lODE 
COMPLEX’16W.WNEWK)FW.DFOFW.POFWPE,WPEW,EPS1L 
COMPl£X’16  FPLUSJ3UMB.WMEW.FOFWMEWCOkR,WSTART 
’  DEFINE  THE  VARIABLE  FOR  THE  OUTPITT  COMMON  BLOCK. 

PARAMETER  (N-iOOO) 

REAL’S  Q2(N).W2(N)jrOOR(N) 

COMPLEX’16  Z(NLP1(N) 

LOGICALINSTAK(N) 

COMMON  /OUTPUT/  PlZQ2,W2JCOOR,INSTAK 
•PRELIMINARIES. 

P1-4JD0’DATAN(IJXD 

TWOPl-lDO’P! 

’  SET  UPTHE  AMBIENT  PRESSURE  LOOP - - 

PMIN-ZOOOS 

PMAX-8.0DS 

PI^-DSQRIXPMIN) 

PMAX-DS(^nXPMAX) 

NPAMBS.16 

’  SET  UP  THE  MODE  LOOP. _ 

DO  IS  MODE^ 

IF(MODEEQ.2)THEN 
TMIN-603.D0 
n  «a? 

TMIN-403XI0 

ENDIF 

DO  17NPRESS-14.NPAMBS 

WCORR  -  (31S1.633184D0,<<7.Sl7S943D0)/(3l52.92m7D0, -86.73400) 
PLOIN-O 
QUAL-30JX) 

PAMB  -  PMIN  (PMAX-PMIN)  *  DFLOAT(NPRESS)/  DF1jOAT(NPAMBS) 
PAMB-PAMB’PAMB 

•HWONTHETEMPERATURECALCULATIONS-- 
DO  10  ITRCH.INT(TMINX1593J0 
EPSIL.IJM 
NTIMES-O 

TRCH-OFLOAirnRCH) 

CALL  CHANGEOaS.TRCH) 

CALLCHANGEOSI.TRGH) 

CALL  CHANCEa56.TRCH) 

CALL  CHANGEOSP.TRGH) 

CALL  CHANCEaS4.TRCH) 

PLOm-PLOTN^l 

CALLTAE(FREMINJFREMAXJIUMFREP0EST) 

POEST  -  POEST’DFLOATIMODE) 
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WRITE  (Mil)  POECT 

1 1 1  PORMATC  RESONANT  FREQUENCY  ESTIMATE  FROM  OPL  •.F9.3) 

•  1F(PL0TN.NE.1)  FREMIN  -  RESFRE‘(l.D0-0^/QUAL) 

•  STARTTHE  COMPLEX  OGENFREQUENCY  ALGORITHM  WITH  THE  DRIVEN  SYSTEM  Q. 

IF  ((QUAL.LT^.DO)JVND.((QUALGT.O.DO)^ND.(PLOTN.LT.7)))  THEN 
DO  100  IFREQ«l,NUMFRE 
IF(NUMFRE^.1)THEN 
FREQ-FREMIN 
ELSE 

FREQ  -  FREMIN  ♦  (FREMAX-FREMIN)»DFLOAT(IFREQVDFLOAT(NUMFRE- 1) 

ENDIF 

FREQ-FREQ*DFIjOAT(MODE) 

W-TWOPI*FREQ 
CALLTAE1(W.P0FW.FPLUSJ»AMB) 

AMPLIT(IFREQ)-CDABS(P1(1)) 

100  FREQUOFREQ)-  FREQ 

•  GETTHE  Q  AND  RESONANT  FREQUENCY  USING  THE  CONSTANT  AMPLIT  DRIVER 
•RESPONSE. 

CALLQUAFAC(AMPLITJREQU.QUALJIESFRE) 

WRITE  (29.1010)  RESFRE.QUAL.TRGH-293.D0 
WRITE  (27,*)  TRCH-293.D0,1JXVQUAL 
WRITE  (28.*)  TRCH-293 J»JIESFRE 
WRITE  (MOlO)  RESFRE,QUAL,TRGH-293.D0 
1010  PORMATC  ROfTE;  RES  FREQ-\F9.3.'  Q-’.P9.3;  DELTAT-’.F9J.’  K*) 

ELSE  iUNEARLYINTERFOUtETOGETASTARTPORQUALANDRESFRE. 

•  U3UALANDRESFRE  ARE  ASSUMED  TO  BEUNEAR  IN  DELTA  T. 

QUAUlJXVQOLIMQOLrM30LD2)»aRGH-TDLDV((TOLDMOLD)‘QOLD2»QOLD) 
QUAL- IDO/QUAL 

RESFRE-  RESOLD  ♦  (RSOU»RESOLD)*fTRGH-TOLDy(TDLD2-TOLD) 

END  IF !  THIS  ENDS  THE  INITIAL  SEARCH  FOR  ASTART  Q  AND  RESFRE 

•  NOW  GETTHE  COMPLEX  EICENFRBQUENCY  FOR  COMPARISON. 

•INITIAL  GUESS  PORTHENEWTONSTECHNIQUECX^ROOTDETERMINATim 

W  .  TWOPI  •  RESFRE  •  (1 JX  •  (OiXOSDOyQUAL} 

WSTART-W 

IF  (((MODE£QJ)jea).(PLQrmj<E.l)).OR.(MODE£Q.l))  W  .  W  •  WCORR 
WRITE  (*,•)  DREAL(WyiWOI»L-DREAL(Wy(2J)0*DIMAC(W)) 
no  CALLTAEIIWJOFWJTUWAMB) 

NTIMES-NTIMES4>1 
WPEW  -  W  •  (1.D0  ♦  EPSIL) 

WMEW  -  W  •  (1  J»  -  EPSl) 

CALLTAEl(WPEW40rWPEJ)UMB.PAMB) 
CALLTAE1(WMEW.P0FWMEJJUMB4»AMB) 

DPOFW-(POFWPE-POFWME)/  (2D0*EPSIL) 

WNEW  -  W  -  KXLDO  •  POFW  /  DPOFW 
EPSn.  -  (WNEW  -  W)«2J>4/(WNEW+W) 

•  TEST-CDABS((WNEW.Wy(WNEW*W)) 

TESTF  -  CDABS(K3FW/FPLUS) 

•  IP  ((TEST  .OT.  1J>4)  .OR.  (TESTF  .CT.liM))  THEN 
IF  (TESTF  .GT.  1J>S)  THEN 
W.WNEW 
GOTO  no 

ELSE  I  COMPLEX  EICB^FREQUENCY  HAS  BEEN  FOUND. 

W-(W  +  WNEW)/2JD0 
WCORR -W/WSTART 
WRITE  (••)NTIMES 
WRITE  (•.•)POFW 

ENDIF  !TEST.CT.ORTESTF.CrCONDrnONAL 

•  GETTHE  Q  AND  RESONANT  FREQUENCY. 

RESFRE  -  DREALCW)  /  TWOPf 
QUAL  - -DREAL(Wy(2X>0  •  DIMAC(W)) 

WRITE  (39,101 1)  RESFRE.QUAL,TRCH-293.D0.PAMB 
WRITE  (37,*)  TRCH-293.DO,l.DOA3UAL 
WRITE  (38.*)  TRCH-293.D0.RESFRB 
WRIIE  (MOll)  RESFRE.QUAl,rRCH-293.D0J*AM8 
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lOll  FORMATC  EIG;  Fff.R.S;  Q\F9.3.'  DT,F9.3.'  5C\'  PAMB-.F10.2) 

•  DETERMINE  IF  ONSET  HAS  BEEN  ACHIEVED;  IF  SO.  GET  DT  AND  RESFRE,  NEXT? 

IFIPLOTN.EQ.  DTHEN 
RESOLD -RESFRE 
QOIJD  -QUAL 
TOLD  >TRGH 
ELSE 

IF  (QOLD^^UAL  iT.  OXW)  THEN 

•  INTERPOLATE  TO  GCTTHE  ONSET  DELTA  T. 

TONSET -TOLD  ♦  CTRGH  -  TOLDVdiX)  -  QOLDQUAL) 

WRITE  (44,*)  PAMB,TONSET-293.DO.MODE 
WRITE  (M 1 17)  PAMB.TONSET-293.DOJl<ODE 

1 1 17  FORMATC  AMBIENT  PRES  •,F10.2.*  ONSET  DELTAT.F9J.-  MODEM3) 

Goroi7 

ELSE  !  RESET  OLD  PARAMETERS  TO  NEW  PARAMETERS. 

RSOLD2.  RESOLD 
QOLCQ  -QOLO 
T0LD2  -TOLD 
RESOLD -RESFRE 
QOLD  -QUAL 
TOLD  -TRCH 

END  IF  !  ONSET  CONWnONAL 

ENDIF  IPLOTNEQ.! 

10  CONTINUE  !  TEMPERATURE  LOOP  INCREMENT 

TNOLD2.TNOLD 
TNOLD-TONSET 

17  CONTINUE  'PRESSURE  LOOP  INCREMENT 

18  CONITNUE  !  MODE  LOOP  INCREMENT 
END 

MAINLINE  USED  TO  COMPUTE  THE  RESONANT  FREQUENCY  AND 
QUALITY  FACTOR  CURVES. 


•  PROGRAM ;  EVALUATE  TAE  PQK  A  RANGE  OF  PARAMETERS. 

•  VERSION  EXPUCTTELYPOR  THE  UMTAE 

•  PAT  ARNOTT.  22  MARCH  1991,  MOD  IM5-91. 

•  DETERMINES  THE  STABUnr  CURVE  FOR  THE  FIRST  TWO  MODES  AS  A  FUNCnON 

•  OF  THE  AMBIENT  PRESSURE.  INCLUDES  FINITE  WALL  POROSITY  IN  THE  STACK. 

•  PROGRAM  USEDTO  DETERMINETHE  Q  AND  RESFRE  OFTHE  UMTAE  VS  TEMPERATURE 


PROGRAM  UMTAE2 

•  VARIABLES  USEDTO  CETTHE  Q. 

REAL’S  AMPL1T(2000)JREQU(2000).QOAURESFREJ1ESOLD.QOLD 
REAL’S  QOLD2,TOLD2J1SOLD2,TOLO 

•  VARIABLES  RETURNED  FROM  TAE 

REAL’S  FREMINJIIEMAXFOEST 
INTEGER  NUMFRE 

’  LOCAL  VARIABLES  TO  THE  MAINLINE. - 

REAL’S  PAMBJ>MIN,PMAX.TMIN 

REAL’S  TRCH.TtWSEr jn.TWOPLTEST,TESrP.TNOLD.TNOLXa 
INTEGER  nRGHJ%OrTNJFREQjmMES,NPAMBSJ<IFRESSJHODE 
COMPLEX’16  W.WNEWJOFWJ>POFWJOFWPE.WPEWEPSIL 
COMPLEX’16  FPLUS4XMB.WMEWJOFWME.WCORR.WSTART 
’  DBTNE  THE  VARIABLE  FOR  THE  OUTPUT  COMMON  BLOCK. 
PARAMETER  (N-3000) 

REAL’S  Q2(N).W2(N),2COOR(N) 

COMPLEX’16  Z(N)J>l(N) 

LOCICALIHSTAX(N) 

COMMON  /OUTTVT/  Pi;Z.Q2,W2aCOOR4NSTAK 
’PRELIMINARIES. 

PI-4I»’DATANaiX)) 

TWOPI-IDO’PI 

’  SET  UPTHE  AMBIENT  PRESSURE  LOOP. _ 

PMIN  -  1.73D5 
PMAX-1.73D5 
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TMIN  -  293.130 
NPAMBS.1 

•  SET  UP  THE  MODE  LOOP . 

DO  18  MODE>2J 

DO  17  NPRESS-1.NPAMB3 

WCORR  .  (3151.633184  D0,-«7.817«943D0V(315Z921237D0.-66.734D0) 

PIjOTTN-O 

QUAU30.D0 

PAMB  -  PMIN  +  (PMAX-PNtIN)  •  DFLOaT(NPRESS)  /  DFLOAT(NPAMBS) 

•  HOP  ON  THE  TEMPERATURE  CALCUUTIONS _ 

DO  10  rTRGH-lNT(TMIN).1393^0 

EPS1L-1JV8 

NTIMES-0 

TRCH-DFLOATHTRCH) 

CALL  CHANGEtt28.TRCH) 

CALLCHANCEa31.TRCH) 

CALL  CHANGEaS6.TRGH) 

CALL  CHANGEa59.TRGH) 

CALL  CHANGEa84.TRCH) 
pijorrN-PijCfrN-..i 

CALLTAE(FREMINJTlEMAXNUMFREJOESr) 

PDEST  -  POEST»DFLOAT(MODE) 

WRITE  (MI  I)  WEST 

111  PORMATC  RESONANT  FREQUENCY  ESTIMATE  FROM  OPL’,F9J) 

•  IF(PIjOTNJ4E.1)FREMIN-RESFRE*(1-DO-OJDO/QUAL) 

•  STAimHE  COMPLEX  E!G»ffREOUENCY  ALGORITHM  WITH  THE  DRIVEN  .SYSTEM  Q. 

IF  ((QUALXTaOO.DO)>ND.(QUAI>CT.OJXJ))  FriEN 
DO  100 IFRBQ-1 JWMFRE 
IF(NUMFREEQ.1)THEN 
FREQ-FREMIN 
ELSE 

FREQ  -  FREMIN  >  (FREMAX-FREMIN)*DilJOAT(IFREQyOFljOAT(NUMFR£-l) 

ENDIF 

FREQ-FREOnSFLOAHMOOE) 

W-TWOin*FRBQ 

CALLTAE1(W^FW4TLUSJ>AMB) 

AMFUTaFREQ)-CI}ABS(Pl(l)) 

100  FREQUOFREQ)- FREQ 

•GETTHEQAND  RESONANT  FREQUENCYUSING  THE  CONSTANT  AMPLIT  DRIVER 
•RESPONSE. 


CALL  QUAFAC(AMPLITJR£QU.QUAU1ESFRE) 

WRiTe  (29,1010)  RESFRE.QUAl,TRCH-293.D0 
WRire  (27,*)  TRCH-293.D0.1  JXVQUAUMODE 
WRHE  (28,*)  TRCH-293J30JIESFREJ«)DE 
WRHE  (•.  1010)  RESFRE.QUAI,TRGH-293.D0 
1010  PORMATf  ROTE;  RES  FREO--.P9.3.’  Q--.P9.3.’  OELTAT-’.F9.3,-  KT 
ELSE  lUNEARLYINTERPOLATETOCErASTARTFORQUALANDRESFRE. 

•  IQUAliANDRESFRE  ARE  ASSUMEDTO  BE  LINEAR  INDELTA  T. 

QUAll-lJX>QOLD*<QOIIM30US)*(TRCH-TOU>)/(aOLD2-TOU;>)*QOLD2*QOLO) 
QUAli- IDO/QUAL 

RESI^  RESOLD  ♦  (ltSOLD2-RESCLO)^(TRGH>TOLDy(rOLIS-TOlD) 
end  !  THIS  ENDS  THE  INITIAL  SEARCH  FOR  A  STARTQ  AND  RESFRE 

•  NOWdETTHECOMPLBCEXIENFREQUENCYPORCOMPARlSON. 
•INnLMiCUESSPDRTHENEWTONSTECHNIQUEOFROOIDETERMINATION. 

IF((MODEEQ.2)JU4D4(QUALQT.100D0).OR.(QUALLT.0.D0))) 

*  TH^  !  THE  ROUTINE  IS  OK  FOR  THE  HIGHER  M(X3E..  CO  AHEAD 
W  .  T^PI  •  RESFRE  •  (1.D0  -  (0D0,0.5D0VQUAL) 

WSTART-W 

IF  (((NtoDEEQJ)uUflD.(PlJOTNNEl)).OR.(MODEJEQ.l))  W  -  W  •  WCORR 
WRITE  (•,*)  DREAL(WVTWOPL-DREAL(Wy(2.DO*DIMAC(W)) 

110  CAliLTAEl(W,FOFWJPLUSJ?AMB) 

NTIMtS.NnMES4.1 

Wl>EVlr-W*(l.D0*EPSIL) 

WMEW  -  W  •  (li»  •  EPSIL) 
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CALL  TAE1(WPEW.F0FWPE.DUMB,PAMB) 

CALL  TAElfVMEW.FOFWWaDUMB.PAMB) 

DFOFW  -  (FOFWPE  -  FOFWME)  /  (IDO  •  EPSIL) 

,  WNEW»W-100.D0»FOFW /DFOFW 
EPSIL  •  (WNEW  -  W)'-ZI>4/(WNEW+W) 

•  TEST-CDABS((WNEW-WV(WNEW+W)) 

TESTF  «  CDABS(FOFW/FPLUS) 

•  IF  ((TE5r  .GT.  1  .D4)  OR.  CTESTF  GT.  1  .D-4))  THEN 
IF  (TESTF  .GT.  l.D-5)  THEN 

\V=.WNEW 
GOTO  no 

ELSt  !  COMPLEX  eGENFREQUENCY  HAS  BEEN  POUND. 

W3(W  +  WNEW)/1D0 
WCORR-W/WSTART 
WRn£(V)NTlMES 
wRrrE(V)  POFw 

LND  IF  !  TEST.GT.  OR  TESTTCT  CONDITIONAL 

•  CETTHE  Q  AND  RESONANT  FR£Q«JENCY. 

RESFRE  -  DREAUW)  /  TWOPI 
QUAL  -  -DREAL(WV(2JD0  •  DIMAC<W)) 

WRITE  (39,1011)  RESFREQUALTRCH-293.D0,PAMB 
WRITE  (37,»)  TRCH-293.DO,1.D(VQUALMODE 
WRITE  (38.»)  TRGH-293.DOJIESFREJHODE 
WRITE  (MOll)  RESFRE.QUAUTRGH-293.DO.PAMB 
1011  FORMATC  EIG;  F(r.F9J.'  Q'.F9.3.-  Dr.F9.3.'  K'.'  PAMB-.F10.2) 

•  DETERMINE  IF  ONSET  HAS  BEEN  ACHIEVED:  IF  Sa  GET  OT  AND  RESFRE,  NEXT? 

IF(PLOTN.EQ.l)THEN 
RESOLD.  RESFRE 
QOLD  -QUAL 
TOLD  .TRCH 
n  q: 

•  IF(QOlI>^UALJ.t.ODO)THEN 

•  imERIXX^TETOCETTraONSETDELTAT. 

•  TONSET-TOID.CTRCH.TOLDV(1DO-QOLD«QUAD 

•  WRITE  (44,*)  PAMB,TONSET-293DOA<ODE 

•  WRITE  (M117)  PAMB,TONSET-293.DO,MODE 

•  1 1 17  FORMATC  AMBIENT  PRES  * Jia2,'  ONSET  DELTAT.F9.3,'  MODEJ3) 

•  GOTO  17 

•  ELSE  !  RESET  OLD  PARAMETERS  TO  NEW  PARAMETERS. 

RSOLDS.  RESOLD 

QOLD2  -QOLD 
TOLD2  -TOLD 
RESOLD -RESFRE 
QOLD  -QUAL 
TOLD  -TRCH 

•  ENDIF  1  ONSET  OONDITTONAL 

END  IF  IPLOTNEQ.! 

ENDIF  ICONDmONALONMOOE£QJANDQ>10a 

10  CONTINUE  ITEMPERATURE  LOOP  INCREMENT 

TNOLDO-TNOLD 
TNOLD-TONSET 

17  CONTINUE  !  PRESSURE  LOOP  INCREMENT 

18  CONTNUE  I  MODE  LOOP  INCREMENT 

END 
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acoustic 


hot  TAE  cold 

heat  (or  heat 

a)  exch  stack)  exch 


single  . arbitrary-perimeter  tube  of  a  thermoacoustic  element 


Figure  1.  a)  Generic  arrangement  used  in  thermoacoustic  heat  engines,  b)  An  exposed  view 
of  a  thermoacoustic  element  consisting  of  a  parallel  combination  of  square  capillary  tubes, 
c)  A  single  arbitrary-perimeter  capillary  tube  for  use  in  a  thermoacoutic  element. 
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Figure  2.  a)  Parallel  plate,  b)  circular,  c)  rectangular,  and  d)  equilateral 
triangular  capillary  tube  geometries 
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^  b) 

LARGEAT:  PRIMEMOVER 
To -AT  To  To+AT 


Hguie  3.  Lagnngian  view  of  a  fluid  parcel  in  a  standing  wave  near  a  boundary. 
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Figure  4.  Magnitude  and  phase  of  the  excess  temperature  between  parallel  plates. 
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thermoacoustic 
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FUNDAMENTAL:  116  Hz 


-5.7  cm 
0  cm 

17.2  cm 

34.3  cm 
68.6  cm 


BACKGROUND  LEVEL:  75  dB 


FREQUENCY  (Hz) 


Figure  1 1.  Spectral  peaks  as  a  function  of  distance  from  the  mouth  of  the  prime  mover. 
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resonant  frequency  (Hz) 
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Hguie  15.  a)  QuaUty  factor,  b)  resonant  frequency  fn*  UM  TAE  with  ambient  pressure  173  kPa. 
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